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SUMMARY 


A general simulation program (6SP) involving nonlinear state 
estimation for space vehicle flight navigation systems is developed 
and used as a basis for evaluating the performance of a Space Tug 
navigation system. A complete explanation of the iterative guidance 
mode (IGM) guidance law, derivation of the dynamics, coordinate frames, 
and state estimation routines are given so as to fully clarify the assump 
tions and approximations involved so that simulation results can be 
placed in their proper perspective. 

So as to make the simulation program as useful as possible, a 
complete set of computer acronyms and their definitions as well as 
explanations of the subroutines used in the GSP simulator are included. 

To facilitate input/output, a complete set of compa table numbers, with 
units, are included to aid in data development. Format specifications, 
output data phrase meanings and purposes, and computer card data input 
are clearly spelled out. 

A large number of simulation and analytical studies are used to 
determine the validity of the simulator itself as well as various data 
runs. Included in the studies are (1) covariance initialization, (2) 
initial offset vector parameter study, (3) propagation time vs. accuracy, 
(4) measurement noise parametric study, (5) deterministic vs. filtered 
run terminal error, (6) reduction in computational burden of an on-board 
implementable scheme, and many others. 
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From toe results of these studies-, conclusions and recommendations 
concerning future areas of important practical and theoretical purpose 
are presented. These use as a basis an extension of the general GSP 
simulator developed to date. 
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I. INTRODUCTION 


A. Definition of Space Tug and Statement of Project Task 

The basic Space Tug vehicle will be used to extend the capabilities 
of the Space Shuttle orbit vehicle. The missions of the Shuttle will be 
limited to low earth orbit, severely limiting mission flexibility. The 
Space Tug concept was devised to extend the mission capability to geo- 
synchronous orbit, defined as a 19,323 nautical mile - 0° circular orbit. 
For maximum flexibility the Space Tug will fly in the cargo bay of the 
Space Shuttle as payload into low earth orbit. Departure of the Tug 
from the Shuttle will then take place in an approximate 160 n. mi. x 
160 n. mi. low earth orbit. This Space Tug with its own payload will 
then fly its desired mission and return to the Shuttle for transportation 
back to earth. 

The primary mission of the Space Tug is to deliver and retrieve 
automated payloads from low earth orbit to geo-synchronous orbit and 
return. Other ^r.ission assignments for the Space Tug are as follows: 

1. Astronomy - Depart from a 160 n.mi. 28. 5‘^ orbit and transfer to 
a 39,000 n.mi. 28.5° circular orbit. 

2. Atmospheric Physics - Depart from a 160 n.mi. 28.5° orbit and 
transfer to an escape trajectory. 

3. Earth Observation - Depart from a 100 n.mi. 90° orbit and transfer 
to a 900 n.mi. 90° circular orbiL 

4. Planetary - Depart from 160 n.mi. 28.5° orbit and transfer to 
an escape trajectory. 
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5. Communication .and Navigation - Depart from 160 n.mi. 28.5° 
orbit and transfer to 19,323 ruini. 0° circular orbit (geo- 
synchronous orbit). Depart from 205 n.mi. 103° orbit (sun 
synchronous orbit) [3]. 

The design configuration of the Space Tug assumed for this project 
is that the vehicle be unmanned, ground based control, reusable and 
launched with payload on a single shuttle flight. It is desired that 
the Space Tug be a low cost vehicle. Under this assumption, existing 
hardware should be incorporated into design and a minimal amount of ground 
support tracking should be employed, necessitating a large burden on an 
on-board inertial navigation-guidance system. 

As stated previously, the primary mission of the Space Tug, and the 
one used as a base mission for this project, is the flight of the vehicle 
from low e,trtf^ orbit to geo-synchronous orbit. The exact definition 
of this mission is: Depart from 160 n.mi. 28.5° orbit and transfer to a 

19,323 n.mi. 0° circular orbit [3]. This mission flight is depicted in 
Figure I-l. The guidance scheme for this orbit transfer had previously 
been adopted by NASA to be a modified Iterative Guidance Mode Law (IGM) 
which is optimal in the sense of minimum fuel [1]. To achieve the orbit 
transfer with minimum fuel, the guidance law will employ a three burn-three 
coast maneuver, instead of a long duration burn. IGM is of the modified 
form in the sense that for the Space Tug the terminal orbit end conditions 
can either be fixed or floating. Another specification of the Space Tug 
flight profile is that the capability of making trajectory changes during 
flight is required. To meet this design specification the guidance and 
navigation scheme needs to have an accurate fix on its position and 
velocity in an inertial coordinate frame at all times. A navigation system 



CIRCULARIZATION AT 
GEO-SYNCH BURN 





Figure I-I. BURNS TO CIRCULAR GEO-SYNCHRONOUS ORBIT 
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with guidance parameters supplied by IGM could not accomplish this task, 
since IGM does not produce the position and velocity vectors of the Space 
Tug, using them instead as inputs. This report describes a simulator 
and navigation system for estimation of the position and velocity state 
vector in an inertial coordinate frame using an extended Kalman filter 
type approach. The overall problem of mid-course guidance and correction 
is involved and for the present purpose it is assumed IGM functions 
properly [2j. 

B. Approach to Problem Solution 

The three position and three velocity coordinates in an inertial 
reference frame make up the 6 x 1 state vector of interest in this 
problem. Because of the nonlinear nature of the vehicle dynamics, which 
include the effects of gravity and aerodynamic drag, the linear Kalman 
filter is not applicable. Therefore, the extended Kalman filter technique, 
which requires linearization of the vehicle dynamics about the nominal 
state estimate, was used. The state estimator scheme uses the linearized 
model in the Kalman Filter to find the best or optimal estimate of the 
state vector at some time, t. The Space Tug's terminal position vector, 

Ry, terminal velocity vector, Vy, and terminal range angle. By, are 
used as inputs to the system. 

The state estimation process uses measurement data supplied from 
various measurement configurations to form a measurement state vector. 

The state estimator then forms a measurement residual vector, which is 
the difference of the actual measurement and an estimated measurement, to 
be used in the nonlinear state estimator to update the state vector. 
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The different measurement, configurations used, for the estimation procedure 
should be kept at a minimum to reduce overall cost, storage, and compu- 
tation time in the on-board computer. However, the measurement systems 
must be capable of supplying enough measurement data to the state estimator 
so that the state estimates meet the error requirements placed on the 
position and velocity vectors. Due to the nature of the Space Tug's 
flight, different measurement configurations during the burn and coast 
phase are considered. 

Three axis accelerometer measurements are made durin§ the burn phase 


of flight. The model of this measurement used in the state estimator 
program is given in (I-l). 



In (1-1) is the sensed acceleration vector, X, is the actual accel- 
oration vector of the vehicle and v is the additive noise vector used 
to model measurement error due to the accelerometer equipment. A 
problem arises as to how the acceleration measurement can be used in the 
Kalman filter. To use a Kalman filtering approach to state estimation 
it is necessary that the measurements taken be a linear combination of 
the states. Since acceleration is not part of the state vector this 
requirement could not be fulfilled in a straightforward manner. It was 
decided that the acceleration measurements could be used to construct a 
measurement vector that is a linear combination of the states. The form 




6 


of the constructed position and velocity measurements are given in 
(1-2) and (1-3). 




X + X. 


At 


At 


(1-2) 


^ ^s * (1-3) 

*■ 

In (1-2) and (1-3) X and X denote the previous position and velocity 
estimates of the state X. More detail as to how the constructed measure- 
ment vector was applied and modified will be given in a later chapter of 
this report. 

Since no measurement data can be derived from the accelerometers 
during the Space Tug's coast phase of flight, a separate measurement 
system is needed. From an investigation of many types of measurement 
devices it was decided to use a horizon sensor and sun seeker during the 
coast phase to obtain the measurements needed by the state estimator. 

The basic geometry for the horizon sensor describing the relation of 
the earth and the Space Tug is shown in Figure 1-2. Using the horizon 
sensor a measurement of the angle (j> is made. Figure 1-3 shows the geometry 
for the sun seeker. The sun seeker is a line of sight device with two 
gimballed axes, from which the angles a and g are supplied as measurement 
data. The measurement vector used during the coast phase of flight is 
given in (1-4). 




SPACE TUG 



Figure 1-2. HORIZON SENSOR GEOMETRY 
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Unlike the measurements during the burn, the measurements during the 
coast can be related directly to the position states of the state vector. 
Since no information concerning the velocity states is given by these 
measurements, the dynamics of the system must be used to reconstruct 
velocity states during the coast. From (1-4) it can be seen that the 
coast measurements are nonlinear, therefore a linearization of the 
measurements must be made. The linearization procedure along with 
analytic expressions relating the measurement angles to the states will 
be given in a later chapter. 

C, General Simulator Package 

In order to demonstrate the capability of the estimation algorithm, 
an entire simulation package of the Space Tug guidance and navigation 
system was developed. By having the entire simulation package available, 
the state estimator can be tested and results verified in detail. The 
general simulator package is a FORTRAN computer program, constructed in 
subroutine form, which simulates all properties of the Space Tug guidance 
and navigation which have an effect on the operation of the state estimator. 
Even though the general simulator package is based on the I6M guidance 
law of the Space Tug, the simulator is general in the sense that guidance 
and navigation systems under the influence of nonlinear dynamics can be 
simulated by the program with only slight modifications. The program 
is constructed so that the initial state vector and the terminal orbit 
conditions are supplied as inputs. These parameters can be changed by 
changing only a few data cards in the program. A detailed explanation of 
the general simulator package will be given in a later chapter. 
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D. Report Outline 

This report is divided into five body chapters with four additional 
supplementary appendices. Chapter II details the development, purpose, 
and operation of the simulator package. Included is a general discussion 
of the IGM guidance law including its input/output characteristics. 

Chapter III introduces the overall nonlinear dynamical problem investigated 
in the Space Tug study and the development of a general navigation state- 
estimation package using the dynamics. Chapter IV introduces and discusses 
the measurement systems for both the burn and the coast, the relationship 
between the measurements and the state vector, and the linearization 
approach used to implement the measurement process in the state estimator. 
Chapter V displays results of simulation studies, along with discussions 
and analyses of the results. Chapter VI covers the findings and con- 
clusions from the work performed, and in addition makes recommendations 
for future work. 

The appendices cover supplementary but important material pertinent 
to a full understanding of the GSP system — Space Tug flight work. 

Appendix A shows the substitutions and partial derivative operations 
required in order to develop the general partial derivative statements 
to be used in <i>(k + 1, k) and H(k + 1). Appendix B covers input/output 
data handling and an explanation of the GSP output. Appendix C 
tabulates overall Tug burn-coast tracking accuracies, initialization 
errors, etc. Appendix D includes a list of symbols and their definitions 
which appear in the GSP computer printout listing. 



II. SIMULATOR DESCRIPTION 


In order to allow for flexible studies with a minimum of prograinming 
burden, a general, nonlinear state estimation type simulator was developed. 
Referred to as GSP (General Stimulation Program), in its present form the 
simulator is a complete description of the gravity and drag forces 
impinging on a space vehicle, as well as the guidance, navigation, 
measurement, and state estimation operations which are involved in space 
vehicle flight. As it is presently implemented there are three 
distinct modes of operation available, (1) deterministic, (2) passenger 
mode, (3) full filter mode. These three modes encompass the major 
conditions of interest in most navigation system studies. For greater 
accuracy, the program is written entirely in double precision. 

The deterministic mode means that the simulator and all operations 
with it (i.e. navigation, guidance, tracking, etc.) are supplied with 
exact position and velocity data. No noise and/or position and velocity 
offsets occur. This condition represents the "ideal case," the best 
that is possible. The second mode, a "passenger" mode, is one in' which 
deterministic data is used to run the simulation as in the first mode, 
but, in addition, a complete state estimator is operated and a 'one step 
propogation' state estimate carried along as a "passenger" for inspection 
of the accuracy of the estimation process. The third, or "full filter" mode, 
uses noisy signals and initially biased position and velocity coordinate 
data to simulate a real-world navigation problem. In this mode a non- 
linear state estimation process is used to provide position and velocity 
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data to the guidance law. The full-filter mode simulates the case 
of an actual flight. 

In the remainder of this chapter details of the overall program, the 
flow of data, input/output data definition, units of data, and external 
user operation of such a simulator package are covered. Available upon 
request from the Department of Electrical Engineering, Auburn University 
is a complete listing of the GSP package. A punched deck utilizing an 
IBM compatabTe 029 punch is also available on request. In Appendix II 
is a complete tabulation of input/output, including an example listing, 
details as to data entry by cards, etc. 

A . General Simulator Operation 

A basic flowchart of the GSP simulator is given in Figure II-l. 

Shown are the major operations the program is designed to automatically 
handle. 

The general flow of data is shown by the arrows in Figure II-l. 
First, an initial state estimate ^(0/0) is input to the program. This 
is the information vjhich is available to the on-board navigation system. 
For the deterministic and passenger mode runs this estimate is the 
actual initial position and velocity coordinates of the vehicle. But, 
for the full filter mode, X(0/0) represents an error value of initial 
state, i.e. 

X(0/0) = XjO) + e(0) 

where X(0) is the actual initial vehicle coordinates and e(0) is the 
initial misalignment error. This error in practice is due to tracking 
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Figure II-l. NAVIGATION SIMULATION FLOW CHART. 
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modeling inaccuracies of the ground-based tracking network. Bounds on 
the allowed size of the elements of c are given in Appendix III. 

This initial state is then input to the targeting program calculations 
along with the rendezvous orbital parameters to determine whether a coast 
or burn phase should occur first. The targeting program is not actually 
included in 6SP. Instead, only the various initial-final trajectory 
parameters for a multi -burn mission to get to geosynchronous orbit are 
available to GSP, the data having been assumed to have come from the 
targeting program. 

In any event, either a burn or coast will occur first. If a coast 
is called for, gravity, drag information is used to develop, on-board 
and in real-time, an equivalent set of state transition and observation 
matrices. The state transition matrix $(k+l,i',) is a piece-wise constant 
6x6 matrix which is needed to propagate the state vector forward in time, 
for use in obtaining a cleaned-up, estimated value of the state at a 
number of seconds later after using observation information. During 
the coast the measurements assumed available were (1) a star tracker, 

(2) a sun seeker, (3) a horizon sensor. The star tracker is used to 
help align the vehicle reference platform by providing angular data. The 
sun seeker provides tracking information in the form of elevation and 
azimuthal angles. The horizon sensor physically measures the half angle 
(j) between the earth's horizon and the centerline from the space vehicle 
through the center of the earth. This particular angle information 
actually determines the vehicle's altitude from the earth's center. This 
"radius vector" plus the two sun angles, is sufficient to yield a unique 
set of position coordinates of the vehicle. 
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Using knowledge of the geometry of the various measurement processes, 
equivalent measurement estimates are made using an equation of the form 

Z(k+l/k) = H(k+1) X(k+l/k) 

where the H(k+1) matrix is used to take state information and construct 
measurements Z. Both the it(k+l,k) and H(k+1) matrices are determined by 
numerical evaluation of analytical partial derivatives programmed into 
the computer. The numerical evaluation is effected by using X(k+l/k) values 
to compute the constants entering $(k+l,k) and H(k+1). 

Next, information using $(k+l/k) and H(k+1), along with the measurement 
data, are input to a nonlinear state estimation process, the end result 
being a "best estimate" (in a mean-rsquare sense) of what the state vector 
is some time At after the last available value. This set of calculations 
involves covariance matrix computations and state update. The output of 

the filter is a new state estimate X(k+l/k+l), 

This new estimate is the only information the vehicle has as to what 

its position is and at what velocities it is traveling. This information 
is then used in place of the original X(0/0) as a new "initial condition", 
only this time at some time instant greater than zero and the process re- 
peats. This closed-loop cyclic operation is clearly depicted in Figure 
II-l, where the output of the filter becomes the input to the guidance law 
(which is bypassed during the coast phase). 

At the appropriate time, main engine ignition occurs and the burn 
phase begins. The initial state value at the beginning of the coast is 
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used in the guidance law (a form of IGM) to compute initial burn angle 
information, where the output of the guidance law are two angles, Xp» 

Xy» pitch and yaw angles respectively. These are used to command the 
main engine gimbals to properly direct the thrust vector. 

Similar to the coast phase, a set of matrices 4>{k+l,k), H(k+1) are 
then developed, based on a series of analytical partial derivatives 
programmed into the navigation system and evaluated at the state value. 
<J)(k+l,k) is evaluated at X(k/k) and H{k+1) at X(k+l/k). The <s>(k+l,k) 
matrix is made up of terms from gravity, drag, and the guidance law. 

It is the complexity of the guidance law which greatly adds to the 
computational burden. A flow chart of the IGM guidance law used is 
shown in Figure II-2. Suffice it to say that the volume of calculations 
involved is very great. A later section of the chapter discusses in 
detail the use, operations, and calculations of the IGM guidance law. 

During the burn a different set of measurements are used than 
during the coast. The measurements in this ise are only relative 
measurements of states, namely the vehicle acceleration in the vehicle 
position coordinates. They are only relative measurements because 
integrating acceleration gives only a relative change of velocity, but 
it never says what the actual velocity is. Integrating the acceleration 
twice gives a relative position change, but never says what the position 
actually is. Using these relative position and velocity changes, a 
measurement is "constructed" from past position and velocity data plus 
the relative changes as determined by the acceleration data. In fact, 
the acceleration data includes not only acceleration due to thrusting, but 
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FIGURE II-2. SIMPLIFIED FLOW CHART OF THE IGM GUIDANCE LAW USED WITH 
THE SPACE TUG (Continued) 
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also acceleration due to gravity and drag. This is accomplished by using 
models of gravity and drag, which in actual practice would be carried 
in an cn-board computer. The accelerations are then added together 
and the relative positional and velocity changes determined. The 
reason drag and gravity model s must be used is that on-board accelerometers 
can measure only relative thrusting, no measuring devices are available 
on-board to measure the accelerations due to these external forces. 

During the burn the measurement matrix H(k+1) is assumed constant, 
equal to the identity matrix, Ig. This is used to form a propagated 
estimate of the measurements, and together, the measurement estimates 
and the measurements are used in the state estimator to come up with a 
one step "cleaned-up" estimate of ^(k+l/k+1). This updated, 'optimal' 
estimate is then used by the guidance law and a closed-loop cycle has 
begun. 

For all three operating modes, i.e. deterministic, passenger, and full 
filter, exact state vector information is available from a fourth-order 
Runge Kutta integration routine. This information is then avail abl.e for 
comparison with the estimated data so as to determine tracking accuracy 
and observation sightings during the coast phase. The means by which 
this is performed in GSP is shown schematically in Figure 11-3. The 
upper flow graph represents the on-board type equations used in the 
navigation system, with gravity, drag, and guidance information going 
into <5, the data based on X information. The arrows indicate direction 
of information (either one Way or two way). The lower line of flow 
represents the exact integration carried out. It is this process that is 
not available in the real problem. It shows that after integrating X{k) 
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to get X(k+1), the geometry of the measurement process, a nonlinear func- 
tion of the state vector, is used to determine a measurement vector, 
to which noise is added. This noisy measurement is then used with the 
on-board filter to determine a new best estimate of state. In order 
to allow for flexible operation, a general gaussian, random noise 
generator with ten independent channels is available in GSP. 

B, Description of the IGM Guidance Law 

Shown in Figure I I -2 is a flow chart of the Iterative Guidance 
Mode (IGM) law being considered for use on the Space Tug as supplied 
by Northrop Services, Inc. This section discusses in detail the purpose, 
plan, and operation of this version of IGM and how it fits into the GSP 
simulator. 

The version of IGM shown in Figure II-2 is a modification of that 
used with the Saturn V SII and SIV-B steering systems. It was selected 
over a Cross Product Steering (CPS) law (as used in the Apollo spacecraft 
CSM and LM orbit transfer steering systems) due to its greater flexibility 
and optimality. It is designed to guide the Tug from low-earth orbit 
to geosynch and back again through a series of phasing orbit plane 
changes. Since basic IGM loses effectiveness when burn times are over 
500 seconds, a modification was provided which slightly biases the time- 
to-burn as a function of the difference (Aa) between the semimajor axis 
of the initial parking orbit and the semimajor axis of the phasing orbit. 
This allows for greater accuracy with burn arcs of over 500 seconds. 

The IGM law is an explicit guidance method which differs greatly 
from previous techniques in that nominal mission profile data is not 
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stored on-board. Instead the IGM law is an iterative step-by-step 
closed-loop solution which enforces the desired end conditions at 
end-of-burn. It is based upon the idea of a linear steering law of 
the form 

X = a + bt (II-l) 

X = a burn angle 

a,b constants 

t = running time with respect to (WRT) a fixed reference 

The law makes use of a flat earth assumption, an average gravity force 
magnitude and direction, and small angle approximations. Performance 
losses (i,e. additional fuel burden) are due to the previous factors 
plus the errors in the predictions of the burn arc and time to go (burn 
time). 

In Table II-l is a list of symbols, with appropriate definitions of 
those variables appearing in Figure II-2. Upon entering IGM, the following 
variables must be known in order to complete the computations: F, I^p, m, 

^g» and either A^, Rj, V.^, Yj fon fixed end conditions 

on Pj, Oj, for floating end conditions. In addition, in both cases 
an estimate of the burn-time- to-go is required even though it changes in 
flight. 

A series of time-varying coefficients are computed and a large 
number of coordinate rotations effected to project present conditions 
to various final coordinate frames and use the differences to drive the 
vehicle. The reference coordinate frame used is the Appolo 13 earth- 
centered launch coordinates (ECLC). 
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Table II-l. List of Symbols for the IGM Guidance Law 


SYMBOL 
F, THRUST 

''ex 

Ijp. SISP 
m 

Tg. TG0 

6t 

X, X 



DEFINITION 
Vehicle thrust 
Exit velocity 
Specific impulse 
Mass of vehicle 
Time-to-go 

Computation cycle time 

Instantaneous position and velocity vector, Apollo 
13 coordinates (z-axis downrange, x along negative 
gravity vector at launch site, y completes right 
nand system) 

Instantaneous gravity vector 
Launch azimuth 
Launch latitude 
Desired terminal incl i nation 

Right ascention of desired terminal descending node 

Desired terminal radius 

Desired terminal velocity 

Desired terminal flight-path angle 

Desired terminal semilatus rectum 

Desired terminal eccentricity 

Desired terminal argument of perigee 

Mass flow rate 

Tau (mass/mass flow rate) 
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SYMBOL 


4>Et 




"y 

at 


DEFINITION 

Instantaneous range angle 

Flight path angle 

Desired terminal true anomaly 

Terminal gravity magnitude 

Pitch gimbal engine burn angle 

Yaw gimbal engine burn angle 

Time-to-go correction due to velocity errors 
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A terminal range angle (j>y is found from 


G = 


cos 

sin 0 |^ sin i 
-sin cos i 


= G X 


-1 




F' T 


0 

cos i 
sin i 

cos (j)|^ 
-sin (J)|^ 
0 


R = 

V = 

Y = sin'\(X • X)/RV) 
V 


sin 0 


n 


N 


-cos ej^j sin i 

cos 0 cos i 
n 


sin (}i|^ sin 
cos ij)|^ sin A^ 
cos Ay 


<I>T = <f>Ej + 1/2 cos Y + ^ cos Yj)(l + 


FK6)T 


(II-2) 


-sin <|)^ 
-cos A^ 
sin A^ 


(II-3) 


(II-4) 

(II-5) 

(II-6) 

(II-7) 

(II-8) 


and Y is the flight path angle, referenced to the local horizontal. This 
on-line computed angle is a function of the instantaneous range angle, 
where the average rate of change of the range angle is 


• - 1 /-? /V cos Y . I X 

***av 'r R.r 


(H-9) 
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i „ times T gives an approximation to the additional angle traversed. 

»v g 

The FK6 is a small adjustment value 

1FK61 < 1 (II-IO) 


which is used to improve overall targeting, hence 


♦ (final) - ♦(present) + - tp^^sent* ' <’ ^ 


A users guide to the development and operation of FK6 is available [4]. 
For a detailed description of IGM concepts and derivation, see [5,6,7]. 

An averaged acceleration vector is found by rotating the present acceler- 
ation values into a terminal frame and comparing it to an equivalent 
terminal gravity vector. 


= 1/2 




(II-12) 


where K is a rotation matrix defined by 


K = '{' G 


41 = 


cos (|)^ 0 sin 

0 1 0 

-sin (j).p 0 cos 


(11-13) 


(11-14) 
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and 5* n. ? represent terminal coordinate frame variables. The 
plane is contained in the desired orbit p1ane« with ^ measured from the 
center of the earth and passing through the orbit injection point, n is 
then normal to the orbit plane resulting in 


Jij = fij = 0 

(11-15) 

is shown in Figure I I -4. 

Because of the definitions of n. c and their earth 

orientation. 

•• p •• •• 

= -u/R^ , ng^ “ ° 

(11-16) 

Sin Yy • Hy 0 t 5y ^y COS Yy 

(11-17) 

?T ® Ry > Hy • Cy 0 

(ir- 18 ) 


Using the definitions, an incremental correction factor aT, representing 
a lumped equivalent of burn time corrections over the cycle time 6t of 
the guidance law (it iterates every 6t seconds), is determined. Terminal 
frame velocity errors are computed assuming the averaged acceleration 
value over the remaining flight from (11-12) is accurate. This is effected 
as follows: 

present 
time 

+ a.„„(terminal frame) • (T„ + aT) (11-19) 

g 


frame) 


= V( terminal frame) 


end-of- 

burn 


estimated remaining 
burn time 




Figure I I -4. SPACE TUG TERMINAL COORDINATE FRAME, 
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Usllig (11-19), velocity errors are computed as 

d ^ kj - ?est (end-of-burn) (11-20) 

where (end-of-burn) is the first entry of (11-19). Similar results 
occur for Ari, a?. The velocity error predictions at end-of-burn are 
really only extrapolated estimates of what type of terminal velocity 
errors will occur if the present thrusting history were to continue and 
if the approximations of (1) averaged gravitational acceleration, and 
(2) Tg + aT actually yields the burn time-to-go. It is true that as 
end-of-burn approaches, all of the approximations become very good and 
terminal homing conditions prevail (small angle approximations hold, 
average acceleration values approach instantaneous, errors tend to zero, 
extrapolation equations become accurate, etc.). 

This extrapolated estimate of the velocity errors at end-of-burn 
in the terminal coordinate frame (ideally the velocity errors should be 
zero) is used to determine relative thrust angles xp, xy» P and y 
representing pitch and yaw respectively. These are computed from 

ip = tan'^ (a|/A5) (11-21) 

iy = tan'^ [An/VuP + A?) (11-22) 

These angles are shown in Figure II-5- A number of parameters A , B , C , 

y y y 

°y’ ^y’ S’ S’ S’ functions of T^, ip, i^, 5, 

n. Ij, n-j-, ?, n, h, 5g* are computed. The end result of interest 

here is that using the listed parameters, steering angle parameters Ky, 

S* S’ S computed. Their purpose is to provide a linear steering 
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law form as in (II-l). In theory they should be constants, but due to 
the large number of approximations used in the development of I6M, they 
are updated every cycle (6t seconds) so as to form essentially a series 
of piece-wise linear steering laws of the form 


Xp' " Xp - + K 2 6t 


Xy' = Xy - Kg + 6t 


(11-23) 

(11-24) 


(the signs on are unimportant as they are a function of the particular 
defining variables). The angles are measured in the terminal (?,n»?) 
coordinate frame and as such must be converted to angle information in * 
ECLC so that the gimballed thruster can be properly oriented to drive the 
vehicle in the proper direction. 

This last set of operations is effected by using Xp'» Xy' to 
determine the direction cosines u, v, w, in the c, n» ? frame. 


(11-25) 


- - 


“ — 

u 

■ 

sin Xp' cos Xy’ 

V 

= 

sin Xy' 

w 


Xp' cos Xy' 


Using the K rotation matrix from (11-13), the corresponding direction 
cosines OCX, DCY, DCZ, in ECLC are found from 


OCX 


u 

DCY 


V 

DCZ 


w 


(11-26) 
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It should be noted that since K is an orthogonal matrix, then = K^. 
This greatly simplifies the computer programming burden. The burn angles 
Xp, Xy can then be computed to be 

Xp = tan‘^ (DCX/DCZ) (11-27) 

Xy = -tan'^ [OCY/ ^DCX^ + DCZ^ ] (11-28) 

which is similar to the previous C-n-5 case. These two angles are the 
output which go to the gimballed thruster servos to orient the thrusters 
properly. 

This entire procedure, ending with angles Xp and Xy. is recomputed 
every 5t seconds. The case just described is one in which fixed end 
conditions are assumed and the terminal point is fixed (time-invariant). 

A similar procedure occurs for the case of floating end conditions. Here, 
Pj. e.|., 0 )^ are fixed but Ry, V^, are variable insofar as they are able 
to cause the proper orbit intercept; they must be calculated on-line, 
however, since the remaining calculations require these values. 

Since this mode essentially results in more unknowns than there are 
equations, a sub-iterative loop within the IGM loop itself is set up to 
match parameters as best as possible. The number of iterations is 
controlled by an externally read in data card variable NCYL. Because of 
the discrete-continuous nature of the IGM equations, the resulting partial 
derivatives for use in the navigation-state estimation process are only 
approximations to the true partial derivatives. The iterative cycle with 
the floating end conditions creates a problem because each iteration adds 
incremental values to T^, velocity differences, etc. so that if more than 
one cycle is run through in the sub-loop, the partial derivatives in the 
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state estimator portion of GSP will not correspond to actual parameter 
values. Both sets of data are required and both must be compatible if 
proper estimation is to occur. This problem is alleviated by setting 
NCYL to a value allowing only a single cycle to occur. The shortcoming 
of this approach is that terminal intercept accuracy can be lost to some 
degree. At present, however, this appears to be the only viable solution, 

C. Basic IGM Optimum Trajectory Concept 

In essence, the IGM guidance law is a technique for meeting guidance 
objectives with a minimum expenditure of fuel. All equations are based 
on Newtonian physics and a point mass vehicle in a vacuum. Two different 
simplified two-dimensional derivations of optimum steering laws for a 
minimum fuel system will be briefly presented and discussed with respect 
to the full IGM law in Section B so as to obtain a clearer picture of 
the reasoning behind IGM. 

1. Flat Earth Equations 

For short range flights of approximately constant altitude, the earth 
can be assumed to be flgt and to have a constant gravity vector. A 
simplified model of the 2-D flight problem is given in Figure II-6. 
Referring to this figure the equations of motion are: 


X 

II 

X 

CO 

X^(0) =0 

(11-29) 

^2 = X4 , 

X2(0) = h^, X^Ctp) = 

(11-30) 

Xj = ap(t)cos 0 . 

^3^°^ ^ ''RI* ^3^^F^ “ ''r 2’ 

(H-31) 

X4 = ap(t)sin e-g. 


(11-32) 


and 




36 



0 X 

REFERENCE 


Figure I I -6, CARTESIAN COORDINATE SYSTEM FOR A FLAT EARTH. 
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where ap(t) is the magnitude of the acceleration on the vehicle due to 
engine thrusting. The auxiliary, or costate equations formed from (II-r29) 
through (11-32) are 

- 0 . ~ (11-33) 

^>2 = 0 * (U-34) 

4»3 = - iii-j, (11-35) 

and 

'1'^ “ “ 4 > 2 ’ (11-36) 

These are obtained by forming the Hamiltonian [8]. Taking the partial 
derivative of h and setting it to zero yields 

0 = tan'^ti/^/i/ij) (11-37) 


Equations (IN33) - (11-36) ar^ readily integrated, resulting in 


'l'-j(t) = 0 
<|/3(t) = iPj(O) 

" '(^4(0) ^ 

Substitution of (11-38) into (11-37) yields 
e = tan'^C^ - Cgt) 

where 


(11-38) 


(II-39) 


'i'4(0) 

'l»3(0) 


'l'p(O) 

C = 

*3(0) 


(11-40) 
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Now (H-39) can be written as 

tan 0 = - C^t (H-41) 

and for small angles tan e = e, so by selecting the flight geometry so 
the angle 0 tends toward zero near the end of burn flight, the control 
law is linear as in (II-l). Since gravity is constant, so are and C 2 * 
The resulting control law from (11-41), 


0 - - C 2 t (11-42) 

is similar in form to (11-23) and (11-24) in IGM. The differences come 
from the fact that in IGM the guidance algorithm is updated every 6 t 
seconds, representing the "life span" of the particular constants 
and also because of the discretized nature of the IGM terms. 

2. Round Earth Equations 

In this section a more realistic approximation of the actual 
flight of a space vehicle is studied. It is assumed that the vehicle 
isl flying about a round earth and that an inverse square law gravity 
field is present. The basic flight geometry is shown in Figure II-7. 

The equations governing flight under these conditions can now be written: 


X3, 


X4. 


- + 3p(t) sin 0 , 

Xl 


and 


* 4 “ 


2X,X. 

- y -'- + ap(t) cos 6 , 

^1 ■■ 


(11-43) 

(11-44) 

(11-45) 


(11-46) 
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Figure I 1-7. COORDINATE SYSTEM FOR A ROUND EARTH. 
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where represents the distance of the space vehicle from the earth's 
center, and X 2 is the range angle defined by Figure II-6. As before, 
the costate equations can be determined using calculus of variations 
(COV) principles to be 


and 




? 2,. 2X-,X. ap(t) 

- ♦3(’'4 * ri) - *4 - - 72 “ »>• 

^ ^1 n 


= 0 , 


2X, 


’^1 + '^4 ’ 


i»4 = - i>2 - 'J^3(2X^X^) + Tp^(2 


\pn - ap(t) (ijj- sin e + y- cos 0 ) 


(11-47) 

(11-48) 

(11-49) 

(11-50) 


(11-51) 


where the boundary conditions X-j(O), Xg(0), X3(0), X^(0), X^(t^), X 3 (t^), 
^^(t^)* known. Again taking the partial derivative 
of H WRT 0 yields 


0 = tan 


-1 


^1 'h 


(11-52) 


which can be put in the form 
tan 0 = 


^1 ’^3 




(11-53) 


To obtain a simple expression for 0 in terms of X-j , X^ would be very 
difficult due to the complex forms appearing in (11-47) - (M-51). This 
clearly shows that, just for the simplified 2-D problem, adding the freedom 
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of a spherical earth with an inverse square gravity law as opposed to a 
flat earth, constant gravity model has added greatly to the complexity 
of the optimum solution form. Considering the full 3-D problem with 
general gravity terms compounds the problem even more severely; hence 
the reason for an on-line, iteratively implemented, piece-wise linear 
type guidance law such as IGM. 

D. GSP Subroutine Operations 

The simulation routine GSP is made up of a number of function sub- 
programs and subroutines in conjunction with a large number of COMMON 
blocks for ease of data transfer. A series of routines is used so that 
each can be optimized to perform a certain repetitive chore with a 
minimum of computer cards and hence computer storage and compile load. 

In this section the subroutines are listed, along with the data 
input to each routine, the resulting data output, and the purpose 
behind each routine. It is hoped that this will give the user a better 
understanding of the flexibility of GSP as well as afford a clear enough 
understanding to be able to correctly modify the package without 
interrupting the proper flow of data. 
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MAIN 


input 

reads e, E, n, T , R of each WRT sun 
(for sun coordi Rates) 

X, Thrust, Igp, m, T^, 

^offset 

Atp Atp print increment coast 

^ ^ print increment burn 

A^. ■ launch coordinates 

ICHECK - mode selection 

1 deterministic 

2 state estimator 

3 passenger 

ISW = burn or coast initialization 

NCYL = # cycles through sub-loop if only final orbit (not position) 
is specified 

NFK6 = whether or not the correction factor FK6 is used 

0 means no 

1 means yes 

Initializes (1) state covariance 

(2) noise covariance 

(3) measurement covariance 

(4) X vector 

T COAST - time coast is to end; added to be used as. an ignition 

equation. This idea suggested by B. Williams when routines 
supplied by NASA failed 


output 


none -intermediate step 
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Purpose 

(1) read In start data 

(2) Initialize states, covariance matrices 

(3) Run GSP 

(4) control output print cycle 

NPRINT 

1 nput 

updated simulation variable values 

States, velocities, range angle, accelerations, gravity, time to burn, 
etc . 

output 

prints out values 
Purpose 

Data output routine 

CONIC 

Input 

present flight conditions 
output 

orbit parameters 
purpose 

as open-loop computation of orbit parameter equivalents 

GRAV 

Input - present vehicle position 

output - 3-D gravity value 

purpose - supply gravity values to simulation 

Fischer gravity model 
see NASA TN D-2691 


March 1965 by Fischer 
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I FUN 


input 


position, velocity vector, time 
output 

accelerations due to thrusting and gravity 


purpose 


to compute accelerations on vehicle as integration 
through fractional integration intervals. This is 
(for simplicity) by computing an approximate angle 


Xp(k+1) 


Xp(k+1) 


At 


X^ 


routine steps 

accomplished 

rate. 


and then assuming that between time (k+1) and (k+2) that x_ is 
constant. With 4th order Runge-Kutta intervals of B, , 
all less than the integration step size At, the corresponding new 
burn angles are approximated by 

Xp(k+l+3^.) = Xp(k+1) + Xp(k+1) • B-j 


input 


ARTAN 


sine, cosine of burn angles, integer selector ISW 
output 

normalized angle data 
purpose 

if ISW = 1, angle normalized between 0 and 2it 
if ISW = -1, angle normalized between -tt and v 


VMAG 

input 

three dimensional vector 
output 

magnitude of vector 
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purpose 

compute magnitude of 3-D vector 

VSUB 

input 

two 3-D vectors 
output 

one 3-D vector 
purpos e 

compute difference of two 3-D vectors 

VDOT 

input 

two 3-D vectors 
Output 

one 3-D vector 
purpose 

compute the dot product of two 3-D vectors 

VUNIT 

input 

3-D vector 
output 

3-D vector 
purpose 

normalize a 3-D vector 
^Snew = 4ld/ l^ldl 
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VCROSS 


input 

two 3-D vectors 
output 

one 3-D vector 
purpose 

compute the cross product of two vectors 

I RUN 

input 

positions, velocities, accelerations, time, integratii:>n step size 
output 

updated positions, velocities, time 
purpose 

Exact integration of system dynamics so that (1) true estimation 
error can be determined, (2) observation sightings of celestial 
objects can be determined for the measurement system simulation 


MATMV 

input 

A - a 3x3 matrix 
B - a 3x1 vector 

output 

C - 3x1 vector 
purpose 

C = A • ^ matrix multiplication 
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ATNA 

input 

sin and cos of an angle x 
output 

angle x normalized to be between 0 and 2tt 
purpose 

normalize burn angle to be between 0 and 2ir 

IGM 

input 

old burn angles, gravity, positions, velocities, time, terminal 
or rendevous conditions 

output 

gimballed thruster burn angles Xp» Xy> for pitch and yaw commands 
purpose 

determine optimum burn angles for main engine 

ATNB 

input 

sin and cos of an angle 
output 

normalized burn angle 
purpose 

normalize an angle so it lies between ±tr. 

ASIN 

input 


sin of an angle x 
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output 

angle x normalized to ±j 
purpose 

normalize an angle x» given sin x> to be between 

ACOS 

input 

cos of an angle x 
output 

angle x normalized to ±|- 
purpose 

normalize an angle x» given cos x> to be between ±|- 

YDOT 

input 

state variable dynamical information; system states 
output 

derivative of state variables 
purpose 

compute of y = Ay + B^. Can be used as part of a separate, on 
board routine for on-line integration of the burn accelerometer 
data 

RUNGE 

input 

system order, derivative of state vector, present state vector, 
present time, integrvition interval 

output 

new state vector at updated time, new time 
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purpose 

use in conjunction with YDOT as part of an on'-board routine for 
on-line integration of the burn accelerometer data. 

MMUL3 

input 

matrix A 6x6 
vector B 6x1 

output 

vector C 6x1 

purpose 

matrix multiply a matrix times a vector C = A-B 

MMUL 

input 

matrix X 6x6 
matrix Y 6x6 

output 

matrix Z 6x6 
Z = X-Y 

purpose 

matrix multiply two 6x6 matrices 

MATINV 

input 

matrix, system order 
output 

matrix 

purpose 


matrix inversion 
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MEASUR 

i nput 

vehicle position, measurement observations, noise generator signals 
output 

perfect and noisy measurement observations 
purpose 

to generate accelerometer and/or optical sighting signals of celestial 
bodies and to make the signal $ noisy. 

FILTER 

input 

present system state, noisy measurements, filter partial derivative 
data from IGM, covariance matrix data, mode of operation (Coast or 
Burn) 

output 

updated state covariance, updated state estimate 
purpose 

nonlinear state estimation 

PROPGT 

input 

present system state estimate )((k/k), linearized state matrix, 
gravity 

output 

one step propagated state estimate, ^(k+l/k) 
purpose 

perform the one step state vector propagation 
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SUN 

input 

exact state vector, eccentricity, eccentric anamoly, mean daily 
motion, reference time (time of pericenter passage), semi-major 
axis dimension 

output 

coordinates of sun WRT the earth in a set of coordinates with the 
plane of the ecliptic lying in one of the coordinate planes. 

purpose 

compute sun coordinates for use in sun seeker tracking measurements 


ATTUDE 

input 


at present, none 
output 

three attitude angles 
purpose 

This routine is to represent the attitude control system and inertial 
platform alignment of the Tug. It is properly called in the program 
cycle; at present, the routine is set up for the platform to be 
aligned. 


PARTLl 

input 

^(k/k), X(k+l/k), all data from IGM calculations, mode type (Coast 
or Burn), gravity 

output 

parts of calculations for <i>(k+l ,k) , H(k+T) 
purpose 

To evaluate intermediate expressions which make up the partial derivative 
calculations. There is a PARTLl and PARTL2, and PARTL3, all performing 
calculations. Three separate routines are desired because of Roll Size 
limitations of the IBM 370/155 system at Auburn University. 
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PARTL2 

input 

same as PARTLl 
output 

same as PARTU 
purpose 

same as PARTLl 

PARTL3 

input 

same as PARTLl 
output 

final calculations for 4i{k+l ,k), H(k+1 ) 
purpose 

same as PARTLl 


input 


Output 


output data from PARTLl, PARTL2. SUN. PROPGT, FILTER. MEASUR, RUNGE 
(if used), IGM. IRUN, IFUN, GRAV. CONIC. MAIN. PARTL3 

output 

partial derivative intermediate calculations, 4>(k+l,k), H(k+1), state 
vector estimate, exact state, exact and noisy observations, covariance 
matrix, state and measurement error residuals 

purpose 


make computation results available; also to compare and analyze 
filter accuracy as a function of time. 
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GRAND 

1 nput 

mean value, standard deviation. Integer channel selection number 
between 1 and 10 

output 

A guasslan random (Gaussian RANDo m) number with characteristics 
of each mean and standard deviation for each channel. 

purpose 

Noise generator for measurement process. 

URAND 

input 

integer representing a channel number between 1 and 10 
output 

a single number 
purpose 

Uniform RANDo m distribution function 

BLK DATA 

input 

none 

output 

none 

purpose 

provide initialization values for the random number generators. 
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input 


CONOBS 


matrices A, B, H and their respective dimensions of the form 
X = Ax + Bu 
Z = Hx 


output 

statement as to controllability and/or observability or the lack 
thereof, and a scalar representing the determinant of a matrix 
result based on Kalman's criteria for linear, time-invariant systems 

purpose 

determine controllability and/or observability or the lack thereof 
of a linear (or the equivelant), time-invariant dynamical system in 
state variable form. 


MATRIX 


input 

none 

output 

none 

purpose 

Actually a package of ENTRY statements which perform operations of 
determining eigenvalues of a matrix, multiply, subtraction, inversion, 
normalization, etc. 

CHARD 

input 

An nxn matrix A and an integer print control character 
output 

Eigenvalues of a real matrix 
purpose 

compute the eigenvalues of a real matrix, symmetric or nonsymmetric 
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DINVER 

input 

A square matrix and its dimensions 
output 

a matrix and a scalar representing the determinant of the input matrix 
purpose 

compute the double precision inverse (Double precision INVERse) of 
a matrix 

POLYRF 

input 

root of a polynomial, coefficients of polynomial 
output 

reduced order root, normalized polynomial coefficients 
purpose 

divide out first order factor from polynomial in order to obtain a 
reduced order polynomial 


POLYEV 

input 

root of a polynomial, coefficients of a polynomial 
output 

reduced order polynomial coefficients 
purpose 

used in conjunction with POLYRF and LEMBRT to find all roots of a 
polynomial 
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LEMBRT 

input 

coefficients of a polynomial and the degree of the polynomial 
output 

one root of the remaining reduced-order polynomial 
purpose 

this routine systematically finds the roots of a polynomial one at 
a time using a simple caging scheme based on D'Alembert's Lemma 
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E. GSP Input/Output 

In this section a discussion is given concerning the input/output 
data, the units employed, any constants needed, units of the output 
data, and what the output data means. As the GSP program is now 
structured, all major data is read in from cards or parameters initialized 
from ENTRY statements from commands at the beginning of MAIN, the cycling 
control routine. 

Following is a list of the variables read in, their definition, units 
required (where applicable), and representative values used at the present. 
In Appendix B details of computer card formats and examples of what the 
output looks like are given. The variables read in are: 

ECC - eccentricity of the Earth's orbit about the sun; unitless. Present 
value = .16727 

ECCAML - eccentric anamoly of the Earth; .0083013 

ANM - mean daily motion of the Earth; rad/sec. Present value = 1.9908 x 
10-7 rad/sec. 

I 

TREF - time of Earth's pericenter passage; seconds. Due to Tack of 
definitive data, arbitrarily chosen as 0.0. 

AMAJOR - dimension of semi-major axis of Earth's orbit about the sun; meters 
Present value = 1.4947x1012 meters. 

ISW - program mode switch; 1 means in Coast, 2 means in Burn. The initial 
value is not really important as the program alters it to fit the 
necessary flight conditions. 

NCYL - integer counter which controls the number of sub-loops performed 
in IGM if floating end conditions are used. NCYL = 0 is used, 
indicating one cycle is perfonned, 

NFK6 - integer mode switch; if 0 it means FK6 = 0 and if 1 it means the 
quadratic curve fit for FK6 is employed. 

THRUST - vehicle main engine thrust. Units Of pounds. Present value = 
15,000 pounds. 
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SISP - engine-fuel specific impulse. Units of seconds. Present value = 

440 seconds. 

WTSTAG - initial weight af Space Tug at simulation initialization time. 
Units of pounds. Present value = 58,500 pounds. 

TF - freeze time before engine cut-off. Used to control certain I6M 
parameters near end-of-burn. Units of seconds. Present value = 

10 seconds. 

TIME - simulation running time referenced to Space Shuttle Earth launch 
time. Units of seconds. Present value = 41975.59 seconds. 

X - 3-D position vector in ECLC. Units of meters. Present values = 
-87143.37 -241404.6 -6673876.0 

DX - 3-D velocity vector in ECLC. Units of meters/second. Present values 
7730.2 -17.55 -105.50 

XOFFST - initial navigation state vector error. Six elements; first three 
are position errors in meters and second three are velocity errors 
in meters/ second. 

DTC - Coast propagation time increment. Units of seconds. Present value = 
60 seconds. 

DPRTC - Coast print increment. Automatically adjusted to be the nearest 
multiple of DTC possible. Units of seconds. Present values 
60 to 500 seconds. 

DTB - Burn propogation time increment. Units of seconds. Values between 
.5 and 2 seconds. 

DPRTB - Burn print increment. Automatically adjusted to be the nearest 
multiple of DTB possible. Units of seconds. Present values 2 to 
10 seconds. 

AZ - launch azimuth on Earth. Units of degrees. Present value = 83.39 
degrees. 

PHIL - launch attitude. Units of degrees. Present value = 28.6 degrees 

TPHIT - not used at present. Use a big number, i.e. 1.0x10^ 

TCHIT - first (of two) characteristics freeze times in I6M. Units of 
seconds. Present .value = 15 seconds. 

ICHECK - program mode selection switch. If 1 then a deterministic run is 
made; 2 the full nonlinear state estimator is employed; 3 the 
passenger mode is run. 
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XINTL (entry statement) - located near the end of FILTER subroutine. 

Used to initialize X(k+l/k+l) by adding the 
offset vector, XOFFSt. 

COVRNC (entry statement) - located near the end of FILTER subroutine. 

Used to initialize the state covariance matrix. 

Only diagonal elements are initially defined. 

Since the diagonal entries represent variances, 

2 2 2 2 2 
a-| -> have units of meter and og 

of (meter/second)^. Present values are 
2 2 ,„8 2 2 ,„6 
°^1 ^3 “ ^ ^5 ~ )U , 

QMATX (entry statement) - located near the end of FILTER subroutine. Used 

to initialize the process noise covariance matri 
Q. Only diagonal elements are defined. Units 
2 2 

of meter and (meter/second) . Present values o 
-> = 2500.0. 

RT - terminal radius scalar from the center of the Earth at end-of-burr^. 
Units of meters. Present value— 6,695,842.0 meters. 

VT - terminal velocity scalar at end-of-burn. Units of meters/second. 
Present value = 8360,4 meters/second. 

THT - Yt 5 desired terminal flight path angle. Units of degrees. Present 
value - 1.254 degrees. 

P02 - semi-latus rectum of orbit the burn is to inject Space Tug into. 

Units of meters. Present value = 7,858,094.0 meters. 

E2 - eccentricity of orbit at end-of-burn. Unitless. Present value = 
.17547. 

ARGPER - argument of the perigee (of the orbit at end-of-burn). Units 
of degrees. Present value = 180.689 degrees. 

DINCL - inclination of the desired burn orbit. Units of degrees. Present 
•value = 29.28 degrees. 

DNODE - longitude of the ascending node of the desired burn orbit. Units 
of degrees. Present value = 99.4 degrees. 

TGB - initial estimate of time-to~go in the burn. Units of seconds. Pre- 
sent value =242.3 seconds. 

RANT - range angle for burn orbit. Units of degrees. Present value = 

171 .028 degrees . 
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WDROP - mass flow rate of the fuel during the burn. Results in a decrease 
in Space Tug weight. 


RMATX (entry statement) 


RCOAST (entry statement) 


- located, near the end of FILTER subroutine. 
Initializes the measurement noise covariance 
matrix R. Only diagonal elements are defined. 


During burn, units are 0 , -> 

2 2 2 
0 ^ Oc (meters/second) . 
‘*2 

0-, -> Og^ = 10.0. 


2 

0 ^ meters , 

Present values are 


- located near the end of FILTER subroutine. 
Initializes the measurement noise covariance 
matrix R during the Coast phase. Only diagonal 
elements are defined. Unitless. Present values 

2 2-4 

are a-| =10 


TCOAST - is a replacement method for use as an ignition equation. It is 

the time the coast is to end. Physically this is to be accomplished 
by comparing rangle angle measurements. Units of seconds. 

Present value = 0.0 (start off with a burn). 



III. FLIGHT DYNAMICS OF SPACE TUG 
AND THE NAVIGATION SYSTEM 


In this chapter the basic dynamical equation considerations are 
investigated and the general nonlinear dynamics of the Space Tug in 
powered and unpowered flight are developed. The equations and derivations 
tacitly assume that the vehicle is a point mass. At this point only the 
3 degree-of-freedom translational problem is considered, although an 
inertial platform with provisions for an attitude control system 
package are included in GSP. 

A. Problem Formulation 

The basic dynamics of the Space Tug flight problem arise from 
Newtonian force considerations, i.e. 

F = ma (III-l) 

For the Tug, assuming solar pressure is negligible, gravity and drag are 

the only external force terms acting on the vehicle. The gravity model 
th 

employed was a 5 order polynomial approximation used on the Saturn V 
Apollo flights [9]. The guidance law is a modified version of IGM, 
which was also used in the Apollo program. It is a practical form of a 
minimum fuel optimal control law with internal terminal position and 
velocity rendezvous control. A basic flow chart of the computational 
operation of IGM was given in Figure II-2. 
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Under the previous conditions then, the system dynamics may be 
written in a fixed, dextral , orthogonal, cartesian frame as 


X 


T/m 

DCX + X 

g 

«• 

Y 

= 

T/m 

DCY + Yg 

"''d 

Z 


T/m 

DCZ + L 


m m 



g 

d 


(III-2) 


where 


OCX) 

DCY > engine gimball angle direction cosines 
DCZ) 

OCX = sin Xp cos Xy + sin Xy + cos Xp cos Xy 

DCY = C2^ sin Xp cos Xy + C22 sin Xy + C23 cos Xp cos Xy 

DCZ = C3^ sin x' cos x^ + C32 sin x^ + C33 cos x^ cos x_^ 

Xp = engine pitch burn angle; in radians. 

Xy = engine yaw burn angle; in radians. 

^ij ~ clenients of inverse rotation matirx (K~^) 

T = vehicle thrust, in pounds. 

^a’ ~ accelerations due to gravity in cartesian coordinates 

^ ^ units of ft/sec2 

m = vehicle mass; units of weight/gravity 

For analysis purposes the vehicle thrust is assumed constant, although 
in practice a stochastic thrust variation can be expected. The vehicle 
mass is assumed time-varying with a constant mass flow rate due to fuel 
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expenditure, although, again, in practice a certain amount of random 
variation can be expected. 

In order to develop a state estimator, (III-2) must be expressed 
in terms of state equations. Since (II 1-2) represents 3 second-order 
systems, a sixth order state vector results. Since the direction 
cosines, gravity, and drag terms are nonlinear functions of X, Y, Z, 

X, Y, Z the overall state equations for (III-2) can be written in the 
general form 

^=f(X) (III-3) 

where 

X^ = (X Y Z X? 

f = nonlinear dynamic vector 

A discrete form of extended Kalman filtering is to be used to process 
the on-board data to obtain the state estimates for the guidance law, 
on-board gravity and drag models. For real-time operation this requires 
the use of a set of 1 inear , discrete equations which approximate (III-3), 
at least over a short time interval. 

To obtain the linearized equations, the function vector X^(t + At), 
representing a propogation disturbance, is expanded in a Taylor series 
about the operating point ]^(t). Such an approach was used by Mehra [10] 
for a nonlinear tracking problem for an Earth re-entry vehicle acted upon 
by gravity and drag forces. Results obtained by such an expansion are 

X(t + At) = X(t) + X + X ^ + H.O.T. (III-4) 
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where H.O.T. represents the higher order partial derivative terms. Using 
(k) to represent time (t) and (k + 1) to represent time(t + At), (III-4) 
can be written as 


X(k + 1) = X(k) + i(k)At + At^ + n(k) (III-5) 

2 ! 

where _n{k) is a random noise vector which is to represent a lumping of 
the H.O.T. of (III-4). But 

i(k) = f(X(k)) (III-6) 

and 


X(k) = 


3f 

n 


dx 

3f 

dt 

" 3X 

k 

k 


f(X(k)) 


Using (III-6) and (III-7) with (III-5) results in 


(III-7) 


X(k + 1) = X(k) + [A(X) + I] f(X(k)) At + ii(k) (III-8) 


where 


3f 

> the Jacobian matrix. 


Defining the state transition matrix as 
$(k + 1, k) = A(X) + I 
then (III-8) becomes 


(III-9) 
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X(k + 1) = l(k) + $(k + 1, k)[f(X)At] + n(k) (III-IO) 

Note that f(^)At represents an estimate of the change in the state vector 
X over the time period At since represents the derivative vector of 
as noted in (III-6). 

The actual analytical computation of A(^) is extraordinarily complex 
and will not be given here due to the length. In Appendix A is detailed 
the analysis and substitution approach used to compute the partial 
derivatives in A()^), and an Addendum is available listing all the partial 
derivatives. The form, however, is straightforward and is given below 

0 
0 
0 

[sficx . 1 + + !!d| 

L 3X m 3X axj 

liDCY . 1 !!a + 

L 3X m 8X 3Xj 

Ucz . I ^ ^ ^ 

|_3X m 9X 3Xj 

#(k + 1, k) = (I + A 

At = propogation time 

Similar to the dynamics in (III-3) there is a set of measurement 
observation equations to be considered. These are of the form 

Z(k + 1) = h(X(k +1)) 



(III-ll) 
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where 


Z(k) is the actual measurement at time k 


h^(^) is the nonlinear vector relationship between the state 
vector and the observation. It is assumed there are no 
dynamics involved in this operation. 


Equation (III-ll) can be expanded in a Taylor series and be placed in 
the form 


Z(k + 1) = Z(k) + H(k + 1) 


. X(k + 1) + v(k + 1) 
X(k + 1) 


(III-12) 


where 


H(k + 1) is the linearized jacobian-type matrix of partial 
derivatives 

y^(k + 1) represents the random round-off effect plus general 
noise in the actual measurement process. 


The actual form of H depends on whether the Tug is in the burn or coast 
mode. Details on the measurement process are covered in Chapter IV. 

Using (HI-9) and (III-12), a basic linear state estimation process 
can be formulated [11, 12] as follows: 

Basic dynamics: 

X(k + 1) = f^(X(k)) 

One step propogation: 

X(k + 1/k) = X(k/k) + $(k +1. k) f{i(k/k))At (IU-13) 
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where the refers to estimated information 
Measurement: 

Z(k + 1) = Z(k) + HX(k + 1) + v(k + 1) (III-14) 

Estimated measurement: 

Z(k + 1/k) = Z(k) + HX(k + 1/k) (III-15) 

Measurement residual : 

Z(k + 1/k) = Z(k + 1) - Z{k + 1/k) (III-16) 

State covariance matrix one-step propogation: 

P(k + 1/k) = $(k + 1, k)P{k/k)f^(k + 1, k) + Q(k) (III-17) 

where 

P(k/k) = E|x(k) X(k)^| (IH-18) 

X = X(k) - X(k/k) (III-19) 

Q(k) = E |ji(k) ii(k)^| process noise covariance (III-20) 

^ matrix 

State estimator measurement weighting matrix: 

K(k + 1) = P(k + l/k)H^[HP(k + l/k)H^ + R(k + 1)]'^ (111-21 ) 

where 

R(k + 1) = E |v^(k + 1) jv(k t 1)^1 measurement error (III-22) 

^ covariance matrix 
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Filtered state covariance matrix: 

P(k + 1/k + 1) = [I - K(k + 1)H]P(k + 1/k) (III-23) 

These filter equations are then iterated along with the process dynamics 
so as to obtain an on-line optimal estimate of the states for the guidance 
law and the drag and gravity models. The overall equations of flight 
are similar to those in [13]. 

B. Simplifications of Partial Derivative Computations 

Of considerable interest with the development of GSP is the reduction 
of computation burden imposed upon the implementable navigation system. 

In particular, the large number of intermediate calculations required 
for the partial derivatives in $(k + 1, k) and H(k + 1) suggest some sort 
of sensitivity and/or simulation analysis. Upon original completion of 
the IGM and measurement partial s, a study was undertaken to determine 
which computations could be easily eliminated without noticeable degrada- 
tion in tracking accuracy. First, a number of the final parameter cal- 
culations in IGM such as were set to zero. Results were very good 

in that very little noticeable degradation of tracking accuracy resulted. 

A number of other terms were investigated but most proved to be too 
sensitive when widely varying mission profiles were considered. Partials 
near the end of IGM are the ones which it is desired to eliminate, for 
by eliminating them, many intermediate computations (of which there are 
many) can be avoided. 

A detailed investigation of this area of study is reported in 
Chapter V. The importance of this work is great, for it is only if a 
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simplified set of calculations requiring a minimum amount of storage are 
available and can be processed in "real time" (with respect to guidance 
calculations) that a practical navigation package can be used on-board 
with the Space Tug. 


C. Simulator Flight Modes 

Figure III-l represents the concept of the position tracking errors 
during the burn phase of the flight, and correspond to the actual 
trajectory of the flight and the estimated trajectory of the flight 
respectively, while represents a desired or nominal trajectory from 
deterministic considerations. This nominal trajectory was obtained by 
using perfect data in the simulation model. T^ and T^ are both corrupted 
by noise due to inherent inaccuracies in the measurement system and 
propagation calculations. Also note that while T_ and T. originate at 

a u 

the same point, T^ is initially offset to represent the incomplete 
knowledge of the initial state vector value. 

During the initial phase of the burn, T and T^ tracked the nominal 
trajectory while staying relatively close to one another. T^ and T^ 
remaining close results in a small normed position error, given by 


|X(f)l[ [Xg(t) - Xg(t)]2 + [Yg(t) - Y^(t)]' 


+ [Zg(t) - Zg(t)] 


n 


(in-24) 


t = t' 


where t‘ is any instant of tinse. 

Equation (III-24) is noted to be time variant and is never expected 
to be zero owing to the fluctuating nature of T^ and T^ as the flight 
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progresses in time. The actual state position vector, [X,(t), Y (t), 

d d 

Zg(t)l, and the estimated state position vector, [Xg(t), Yg(t), Zg(t)], 
are both in ECLC. 

Figure III-l also shows the desired, actual, and estimated intercept 
points. A, B, and C respectively, after the burn phase, evaluated at 
t = tg^fj burn* though the normed error is "small", this 

does not guarantee that the actual intercept point, B, is close to the 
desired intercept point, A. However, this error is the only information 
available during an actual flight. For this reason, the flight is broken 
down to three "burn-coast" phases in order to minimize the error between 
the desired end point and the actual end point. 

Since the desired goal is not to follow the desired trajectory 
exactly, but to reach the terminal point with minimum fuel cost, it is 
intuitive that one should try to keep the normed error "small" and there- 
fore keep the terminal error between A and B small. 
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IV. MEASUREMENT SYSTEM 

In this chapter is presented the details and derivations of the 
observation system used at present for the on-board inertial guidance 
system for the Space Tug. The development of the measurement equations 
relating the states to the actual measurements are presented, along with 
the appropriate partial derivative calculations for the development of 
H(k + 1). 

A. Development of Observation Equations 

As in all Kalman filtering schemes, there must be at least one 
measurement taken for implementation in the filter. These measurements, 
while they may not be actual states, should be able to be manipulated 
to yield state information, and they must be causally related to the 
states. How many and what kind of measurements to be taken plays an 
integral part of any state estimation program. 

Actually, the measurements are not used in their raw form but rather 
are used to make up the "measurement residual". In equation form, 

Z=Z-Z (IV-1) 

where Z^ = measurement residual 
Z = actual measurements 
Z^ = predicted measurements 




The development of the predicted measurements, which parallels closely 
that of the actual measurements will be discussed later. 

The measurement system can best be discussed if divided into two 
areas, (1) the burn and (2) the coast phases of the flight. Each 
makes use of different numbers and types of measurements. 

During the burn, strapdown accelerometers are used to obtain the 
actual measurements. These devices measure accelerations along the 
three principal axes of the tug centered coordinate system. The 
governing equations are: 

z(2) = Xy + Yy (IV-2) 

2(3) =• + Y^ 

where X , H , X = acceleration's along the three principal axes of the 
^ ^ tug centered system 

Yx> Yy. Y^ = noise 

Before being used, these measurements must be transformed into Earth 
Centered Launch Coordinates. A description of the different coordinate 
systems and the rotations require^ between them will be given in the 
discussion of the coast phase measurement scheme. 

As is often the case, the states are not measurejd directly. There- 
fore, they must be reconstructed from the available measurements, i.e. 
the measurement must be transformed into the "^ate estimate. During the 
burn, the state estimates are obtained by integrating each' accelerometer 


measurement. 
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* h 


(IV-3) 

(IV-4) 


where ^VAIL ~ available velocity measurement 

)^VAiL “ available position measurement 

• ^ 

-O’ ^ = available initial condition estimates 

Since the integration yields only the incremental change in the 
state value, the initial conditions, ^ and must be added. This in 
itself presents a problem as only the estimates for these quantities 
are available. 

The dependence on previous estimates in the current measurements is 
a totally undesirable situation that for the moment must be tolerated. 

It in fact violates the Kalman filter assumption that the measurements 
must be independent of one another. A further discussion of the problem 
of these indirect measurements will be given at the end of this chapter. 

Because the accelerometers no longer give useful information, the 
coast phase must make use of a completely different measurement scheme. 
Instead oV state measurements, only indirect positional data can be 
obtained using devices to measure angles between the tug and various 
reference points. The basic geometry of the coast measurements is shown 
in Figure IV-1 , where Rgj, and are measured in star tracker 
coordinates [14]. 

The desired vector in the configuration is To find it, however, 

will require some knowledge of all three vectors. 
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Even though all of the measurements are taken in the same coordinates, 
they must be transformed to Earth Centered Launch Coordinates to be used 
in the filter. This is accomplished by a series of coordinate rotations. 

There are four basic rotations needed to transform the information 
in successive stages between the five coordinate systems used. In 
addition, all but one will be time-varying. The first coordinate system 
is the one in which the measurements are actually taken, the star- 
tracker coordinates [15, 16]. The second coordinate system is centered 
at the platform upon which the measurement devices are mounted. As the 
vehicle flies, the measurement devices must be able to move to remain 
sighted on the proper object. Since the platform is in a fixed position, 
this rotation will be time-varying. Next comes a fixed rotation from 
the platform to a tug centered system. The three axes of this system 
coincide with the three principal axes of the tug. This rotation accounts 
for the fact that the axes of the platform system are not aligned with 
those of the tug centered system. From there, the measurements must 
be rotated through the Euler angles, ¥, e, and (j>, to an earth centered 
inertial system. The final rotation puts the measurements into Earth 
Centered Launch Coordinates (ECLC) [17, 18]. In this coordinate system 
X points along the negative gravity vector at the launch size, Z points 
downrange, and Y completes the right-handed, Cartesian set. 

The rotations can be lumped together into the following matrices: 

Packer " ta][3][R][Aji^ (iv-5) 

SYSTEM 
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where a = measurement device to platform rotation 
e = platform to tug centered rotation 

R = tug centered to inertial earth rotation (contains the Euler 
angles) 

A = inertial earth to ECLC rotation 
For the problem at present, the measurement system will consist of 
rigidly mounted devices whose reference axes coincide with those of the 
tug centered system. Thus, the fir^t two rotations are no longer required, 
and a and 3 become identity matrices. 

The rotation to the inertial earth frame is really three separate 
rotations involving the three coordinate axes. For purposes of this 
study, it will be assumed that a 2-3-1 rotation sequence should be used. 

The matrices needed to carry out the rotations are given below. 




Ry = 


<|) 


cos 

0 

0 

-sin 0 

0 


1 

0 

sin 

0 

0 

cos 0 

cos 

¥ 

sin ¥ 

0 

-sin 

¥ 

cos ¥ 

0 

0 


0 

1 


1 0 0 

0 cos <p sin 4> 

0 -sin ^ cos ij) 


(IV-6) 


(IV-7) 


(IV-8) 




where R„ corresponds to rotation 2, R,„ to 3, and R^ to rotation 1. 

X 0 T q) 


y 
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A 2-3-1 rotation requires that the matrices be multiplied in reverse 
order, i.e., 1-3-2. Performing the indicated multiplication yields the 
following matrix 


cos 4' 

-COS <j) sin V cos 9 + sin ({> sin 4- 
sin 4' cos 9 sin 41 + cos 4 > sin 9 


sin 4> 

cos <p cos f 
-sin (j) cos 4” 


-sin 9 cos 4' 

cos (j> sin 4" sin 9 + sin (j) cos 9 I (IV-9) 
-sin 4) sin 4< sin 0 + cos 4> cos 9 

The final rotation to Earth Centered Launch Coordinates is carried 

v» 

out through the elevation angle, 41 ^^, and the azimuth angle, A^, as shown 
in Figure IV-2. The two rotation matrices required are shown below 


(IV-10) 



cos 

4 >l^ sin 



0 

R<i>L = 

-sin 

4 >l^ cos 

<i>L 


0 


0 

0 



1 


- 






P 




_ 


1 

0 


0 



0 

sin A^ 

- 

■cos 



0 

cos A, 


sin 

A, 


- 

z 


z. 


(IV-11) 


Multiplying the matrices in the order shown yields the proper matrix 


denoted by A in (IV-5) 
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os <{ij^ sin (j)j^ sin ’ -sin cos 
in cos *L sin A^ -cos ({.j_ cos A^ (IV-12) 

0 cos A^ sin A^ 

Premultiplying (IV-12) by (IV-9) yields the total rotation matrix, 

> 3 




cos 0 cos <j)|^ 


-sin sin ¥ 


cos ¥ cos 0 sin (|), 


sin + sin ¥^os (j>^^ sin 


-cos A_ sin 0 cos ¥ 


-cos ¥ cos 0 sin <})j^ 

cos A^ - sin ¥ cos i|)j^ 

cos A^ - sin A^ sin 0 cos ¥ 


-cos <|) sin ¥ cos 0 
cos (|)^ + sin (j) sin (j> 
cos <|)j^ - sin ({)|^ cos <p 
cos ¥ 


-COS (}) sin ¥ cos 0 sin (jij^ sin A^ 

+ sin (}), sin A sin <|) sin 0 + 

I. z 

cos <() cos ¥ COS 4|^ sin A^ 

+ cos A^ COS (() sin ¥ sin 0 
+ cos A^ sin (p cos 0 


- sin (jij^ cos A^ cos ()> sin ¥ cos e 

- sin (j)|^ cos A^ sin ({> sin 0 

- cos (j)j^ cos A^ cos (ji cos ¥ 

+ cos (}i sin ¥ sin e sin A^ 

+ sin A^ sin (}. cos 0 


(IV-13) 


cos (|)|^ sin ¥ cos 0 sin sin A^ sin ¥ cos 0 sin t 

sin ♦ t cos cos ♦ * ® \ ^ 

~ sin ()) cos ¥ cos sin A^ 

sin 0 + sin (j). sin (j) 

- sin (j> sin ¥ cos A^ 

+ cos (|) cos 0 cos 


- sin ¥ cos 0 sin p sin (j>|^ cos A^ 

- sin (|)j^ cos A^ cos 4^ sin 0 
+ sin (j) cos ¥ cos i{>j^ cos A^ 

- sin A^ sin 41 sin ¥ sin 0 
+ cos (j> cos 0 sin A^ 


cos ¥ 
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Having defined (IV-5) can be rewritten 

■Cracker " ^clc 

SYSTEM 

From Figure (IV-1), it can be seen that 

Eer ' % ■ hB 

Writing l^j in its major components yields 
hi ^ ^BT^ 


(IV-14) 


(IV-15) 


(IV-16) 


where I, J, K are the unit vectors of the star tracker system. 

These unit vectors are related to the Earth Centered Launch 
Coordinates by the transformation described above. In equation form 



where i, j, k are the unit vectors in ECLC . 

By using the relationship in (IV-17), the elements of can 

further be expanded as 

^BT ^ ^11 ^^ET ■ ^EB^ * *12*-% “ ^EB^ ^ ^13^^ET " ^EB^ (IV-18) 


Ybt = *21^^ET ■ ^EB^ ' ^EB^ ^23^% ' ^EB^ (IV-19) 

% = Slf^ET ■ ^EB^ " *32[Y^, - Y^g] + tggiZg, - Z^g] (IV-20) 

where the are elements of the transformation matrix T^ 

Xet» Ygy, Zgy are the components of ^ 

Xeb» Ygg, Zgg are the components of ^g 
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The calculation of Y^y, and Z^j, the earth to sun coordinates, 
will be shown later after a discussion of the measurement hardware. 

The three measurement devices used during the coast are the horizon 
sensor, the star tracker, and the sun seeker. Each device contributes 
at least one non-redundant angle measurement to be used to compute a 
position fix. 

The horizon sensor measurement is shown in Figure 1-2. 

Defining = esc <p yields 

Z3 = R/r^ (IV-21) 

where r„ = radius of the earth 
e 



Since r^ is known, by measuring 4 , R can be computed. Thus the 
horizon sensor is used mainly to obtain a radius fix. 

The sun seeker and star tracker provide essentially the same angle 
measurements. The only difference is the actual object sighted upon, 
either the sun or a star. The data is manipulated in the same way, 
i.e., the same transformations are used. Only the sun seeker will be 
shown here. 

The sun seeker provides the additional two angles, elevation and 
azimuth, necessary to reconstruct the position vector in the filter. 
Together with the angle measurement from the horizon sensor, they make 
up the measurement vector. The angles that are measured are shown in 
Figure 1-3, where a is the azimuth angle and g is the elevation angle. 

Letting Z, = cos a 

(IV-22) 

Z 2 = cos e 
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yields 

(IV-23) 

(IV-24) 

In summary, the measurement vector during the burn is 

Z(k + 1) = [X Y Z M 2]’’^ (IV-25) 

and during the coast it is 

Z(k + 1) = [cos a cos s CSC <}>]^ (IV-26) 

It should be remembered that the measurement vector is not composed 
of the actual hardware measurements but, rather, quantities derived from 
them in a manner previously discussed. To review, the actual measure- 
ments during the; burn are accelerations while during the coast they 
are the angles themselves and not the trigonometric functions of them. 

Returning now to the previous discussion, the components of the 
earth to sun vector, are generated by using Kepler's equation 

n(t - T) = E - e sin E (IV-27) 

where n = mean daily motion of the earth 

E = eccentric anomaly ' ^ 

e = eccentricity 




T = time of peri center passage 
t = running time variable 


An examination of equations (IV-18)-(IV-20) will show that the 
measurements are now in terms of the desired states, X^g, Y^g, and Z^g, 
with only those states as unknowns. This is the desired configuration 
for use in the filter. 

An expression for the predicted measurement, which was mentioned 
at the outset of this section, can now be obtained. This predicted 
measurement is, as one would think, the value of the measurement at time 
k + 1 based on the predicted state at time k + 1 given data through time 
k. In equation form. 


Z(k + 1/k) = H(k + l)i(k + 1/k) + E[V(k + 1)] (IV-28) 

where Z^(k + 1/k) = predicted measurement 

X(k + 1/k) = predicted state at time k + 1 based on data through 
time k 

E[V_(k + 1)] = expected value of the measurement noise (assumed 
to be 0) 

B. DEVELOPMENT OF H(k + 1) 

In this section the measurement equation approach of the last 
section is used to develop the measurement matrix, H, which is required 
in the state estimation process. The detailed partial derivative 
computations are fully explained so that a clear understanding of the 
approach utilized may be obtained. 

H(k + 1) is the matrix that relates the measurements to the states, 
i.e., it transforms the state vector into a set of measurements. It is 
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of dimension m x n where m is the number of measurements and n is the 
number of states. 

The derivation of the H(k + 1) matrix can now be shown in light of 
the recent insight into the measurement system's operation. During the 
burn phase, the calculation of H(k +1) is a simple matter. Since the 
elements of the measurement vector are states, the predicted measurement 
vector must also be the predicted states. Therefore, Z(k + 1/k) equals 
X(k + 1/k) and from equation (IV-28), H(k + 1) must be the identity matrix 


H(k + 1) = 


1 0 0 0 0 0 
0 1 0 0 0 0 
0 0 1 0 0 0 
0 0 0 1 0 0 
0 0 0 0 1 0 
0 0 0 0 0 1 


(IV-29) 

cV 


During the coast, however, the calculation of H(k +1) is not so 
direct. The measurement vector in this case consists ni the cosines or 
cosecants of the angles sighted by the horizon sensor and sun seeker (star 
tracker). Some means must be found to transfer the predicted states into 
the set of predicted measurements. 

The first step is to recall the fact that given equations (IV-18)- 
(IV-20) and using Kepler's equation to generate X^j, and the 
measurements are now in terms of the states with only those states as 
unknowns. One approximation for the predicted measurements can be obtained 
by expanding the expression for Z(k + 1/k) in a Taylor series yielding 
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8i(k) 

Z(k+l/k) = Z(k) + -^ 


(X(k+l/k)-X(k/k))+H.O.T. 

X = X{k+1/k) 


(IV-30) 


where H.O.T. are the second and higher order terms hereafter neglected 
as insignificant 

Z(k) = predicted measurement at time k 


Thus, the predicted measurement will be of the form 


^NEW -OLD ^ - 


(IV-31) 


Relating (IV-28) to (IV-30), it can be seen that 


H(k + 1) = 


8Z(k) 


9X 


(IV-32) 


Performing the indicated partial derivatives on (IV-21), (IV-23), and 
(IV-24) yields 


^^1 _ ry 2 Y 2 

3X ■ ~W ^BT ^ 


(x 




') (’'bt' 


(IV-33) 


9Z-1 7 7 

[XoT + Y 2] 


3Y 3Y "''BT 'BT 


r(^ 


3X 


BT 


3Ypt\ / 2 

^BtI^BT ^BT ^T~J y^BT ^ ^BT ) 

(IV-34) 


3Zi 3X 


BT 2 


»r = TT »bt‘ ^ ''b/J’’* - *bt(xbt ^ ^ ('‘bt' " 




(IV-35) 
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fl= 

si 

3Z-I 3Z-I 

_A=^= 0 

3Y 3Z 


1 

<m| 

3X|3t 

^hl 3X ^BT 

3X ^ 

3X 

4 /x ^ + Y ^ 
^BT 

^X^+Y ^+Z ^ 
^^BT ^BT % 


SXqT ^^RT 2 2\"^ 

^^BT 3X ^BT 3X ^BT 3X ^^^BT ^ 


2 2 2 
^^BT ^BT ^BT ^ 

3Zg 

3Xrt 

rv “■ + Y 

^''BT 3Y ''BT 

^^BT, 
3Y J 

3Y 

^ ^b/ 

4 /x ^ + Y ^ + Z ^ 
y^BT ^BT ^BT 


(% 3?""''bT 

8YnT 2 2 \^ 

Dl.y d 1\ /y ^ + Y 

3Y "^BT 3Y ^ ^^BT ^BT ' 

- 


(Xb/ " ''bt^ * 

^BT ^ 

3Zg 

3Xr,T ^^RT 

fX + Y 

3Z ^ 3Z ^ 

3Z 

Jv 2 + Y 2 

>''BT ^BT 

■ 4 /x ^ + Y ^ + Z ^ 
Y BT ^BT >BT 


3XrjT ^2dt 0 2 

^^BT 31 ^BT 3Z ^BT 3Z ^ ^ 


l^^BT ^ ^BT * hi ^ 

3Zg _ 3Zg _ aZg _ ^ 

3)J 3Y 3Z 


(IV-36) 


(IV-37) 


(IV-38) 
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^^3 _ 1 ^BT 
3X " r 


V 


e ^/v ^ + Y ^ + Z ^ 
^BT '^BT ■ ^BT 


^ _ J 
aY " r 


'BT 


V 


e ^ + Y ^ + Z ^ 
^BT ^BT ^BT 


^^3 _ 1 

aZ r 


-BT 


{■ 


e V ^+Y ^+Z ^ 

^BT ^BT ^BT 


ax a"? 3^ 


Putting (IV-33) over a common denominator 


aZ 

aT 


(XpT^ Y 2) . x,^(X, 


+ Yot -^) 


1 = ax ^^BT 'BT ^ ^BT^^BT aX 'BT aX 

(X ^ + Y 
^^BT ^BT > 


V 2 ^^BT . V, 2 „ 2 ^^BT 

AnT «. \/ * » DT »%. V “ Af 


3Y 


- X.^Y. 


BT 


_ "BT aX 'BT aX "BT aX "BT'BT aX 


^^BT ■*■ ^BT ^ 


372 


Y 2 ^^BT Y Y 
"^BT aX ■ ^BT’^BT aX 


* ^BT ^ 


3/2 


(IV-39) 


(IV-40) 


(IV-41) 


(IV-42) 


(IV-43) 


{IV-44) 


From (IV-18) and (IV-19) 
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“ax ni 


(IV-45) 


“ax ^21 


(IV-46) 


Substituting (IV-45) and (IV-46) into (IV-44) 


3^1 ^BT ^BT ^21 ' ^BT^ ^11 
aX “ 9 9 3/2 

^BT ^ 


(IV-47) 


Using the same manipulations (IV-34) and (IV-35) can be written 


az 


aY 


1 _ 


7 SXpT 9^ 

yj C D I V V 

''bt aY ■ ■'bt ’bt aY 


BT 




3/2 


(IV-48) 


aZ^ 

aT 


? 9XpT 9Y 

yj C D I Y Y 

^BT az ■ ^BT '^BT aZ 


BT 




3/2 


(IV-49) 


From (IV-18) and (IV-19) 


ax 


BT 


aY 


= -t 


12 


(IV-50) 


aX 


BT ^ 

aZ 


= -t 


13 


(IV-51) 


aY 


BT _ 


aY 


= -t 


22 


(IV-52) 



91 


!!bi_ , 

9Z " "^23 


(IV-53) 


Substituting {IV-50)-(IV-53) into (IV-48) and (IV-49) 


9Z^ _ *22 hj Ygx ' *12 ^BT 
3Y ■ 9 5 3/2 


(IV-54) 


aZ] _ t23 Xg^ Yg^ - t^3 Yg^ 

9Z ; T37I 

^BT ^ 


(IV-55) 


Putting (IV-36) over a common denominator 


9Z„ (X 


9X 


BT 


9Y 


+ Y 


9X 


2 _ '"BT 9X 'BT 9X ^^''BT 'BT ^BT 


^^)(X 2 + Y 2 + z 2) 




3/2 


'BT 


[(Xdt + Yrt ){X 


9X 


BT 


9Y 


+ Y. 


BT 


9Z, 


+ Zr 


Dl 


BT 'BT '^''BT 9X 'BT 9X ' BT 9X 


)1 


^^BT ^BT ^ ^^BT^ ^BT^ ^BT^^ 


3/2 


(IV-56) 


Expanding the numerator and cancelling terms 


y 3 ^^BT y Y 2 ^^BT y ,2 ®^BT . y Y 2 ^^BT 

'^BT 9X ^BT '^BT 9X ^BT ^BT 9X ^'BT ^BT 9X 


+ Y 2 ^^BT Y ■'2 ^^BT y 3 ^^BT y 2 y 

* • r»T «. u ’ 1 iSt ^n*r r.\j ” A|-»t _ v ** T r 


9Y 


BT 


'BT 9X 'BT ^BT 9X “ ''BT 9X BT 'BT 9X 


(IV-57) 


Y 2 7 ^hj 

AriT ^ 


- Y 


^^BT Y 3 ^^BT Y 27 

* I rfcT V “ TfsT 


9Z 


BT 


BT '■BT 9X 'BT ''BT 9X 'BT 9X 'BT ''BT 9X 
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Combining (IV-56) and (IV-57) 


3X 


aZ, [X 


BT 


3X 


2 _ ^''BT 3X 


+ Y — ^IZ ^ - fX ^ + Y ^1Z ^ 

'^BT 9X J^BT ^^BT ^BT ^^BT 3X 


(^BT ^ ^ Yrt ■*■ ^RT ) 


3/2 


BT 


BT 


'BT 


-BT 


From (IV-20) 


^^BT _ . 
3X “ ’^31 


Combining (IV-45), (IV-46), (IV-58), and (IV-59) yields 


3Z2 (-t^^ Xg^ - + {^BT^ * %^^^BT^^31^ 


3X 


^^BT^ ^ ^BT^^ ^’^BT ^ ^BT * ^BT^ 


3/2 


By analagous reasoning (IV-37) and (IV-38) can be written 


3Z, (X 


3X 


BT 


3Y 


2 _ '''BT 3Y 


+ Y ^)Z 2 - (X 2 + V 2)7 ^ 

^BT 3Y ^^BT ^^BT ’^BT ^^BT 3Y 


^'bt' " ^Bt')^' ^^BT^ " Yg/ . Zg,2) 


J/T 


aZ, (X 


3X 


BT 


2 _ '''BT 3Z 


+ Y 


'V)z 2 _/x 2 Y 2w '^BT 

;Z.RT VAgj T TpjY ;4p 


BT 3Z ^"BT 


3Z n o % o o o 

^^BT * ^BT ^ ^^BT ■*■ ^BT * hi ^ 


BT '^BT 3Z 

W2 


From (IV-20) 


3Z 


BT _ 


3Y 


= -t 


32 


(IV-58) 


(IV-59) 


(IV-60) 


(IV-61) 


(IV-62) 


(IV-63) 
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= t 

3Z "^33 


(IV-64) 


Substituting (IV-50)-(IV-53) , (IV-63), (IV-64) into (IV-61) and 


3Y " 9 P Js o o o 3/2 

^^BT ^BT ^ ^^BT ^BT ^BT ^ 


^^2 _ ^'^13^BT'^23^BT^^BT ' ^^BT^ ^ ^BT ^^^33^^BT 
3Z o p % 999 3/2 

^^BT ^BT ^ ^^BT ■*■ ^BT ^BT ^ 


All of the partial derivatives that are necessary to implement 
H(k + 1) in the filter have now been derived. 

Expanding equation (IV-32) 

3Z-| 3 Z 2 ^^3 

3)T" w ~ w 

3Z-| 3 Z 2 ^^3 

^ 9Y” 

3Z^ 9 Z 2 aZg 

3Z“ W W 

An examination of the H(k + 1) matrix, or the partial derivatives, 
will reveal a problem peculiar to the coast phase. All the partial 
derivatives with respect to the velocities are zero, indicating that 
the measurements have yielded no velocity information. In other words, 
no velocity measurements, either direct or indirect, have been taken. 
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Thus, the accuracy of the velocities rely solely on the accuracy of the 
partial derivatives from the dynamics. 

The concept of direct and indirect measurements, which was mentioned 
briefly in the proceeding paragraphs, deserves some consideration. 
Actually the concept is self-explanatory; it is its effect on the state 
estimator that is of interest here. A direct measurement is one in 
which the measurements are directly related to the actual states. A 
single transformation matrix would at most be required to transform the 
data into a usable form for the filter. An indirect measurement means 
that no such direct state information is available. The value of the 
state measurements must be inferred from a combination of the present 
measurements, past measurements, and an initial estimate. An examination 
of equations (IV^3), (IV-4), and (IV-16) shows that, because of the 
initial conditions contained therein, all of the measurements are 
indirect ones. 

The effect of indirect measurement on the normed error can be seen 

without much difficulty. If at the start of a cycle, either burn or 

coast, the initial conditions are in error, then the predicted and actual 

measurement will also be wrong by at least the amount of error in the 

initial conditions. Of course, any errors introduced by the measurement 

system itself will only add to the overall error. Since the initial 

s t 

conditions for the (k + 1) measurements are the measurements at time 
k, each measurement will contain the initial error. Thus, under these 
constraints, the best that a filtering scheme can do is to prevent the 
growth of the error. It appears that it cannot be completely eliminated. 
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C. System Observability 

Having designed the measurement system, there must be an assurance 
that this system will be of use. This assurance comes from considering 
the observability of the system, that is, whether or not the desired 
states can be inferred by observing the output of the system. A 
necessary and sufficient condition for the system to be observable is that 
the n X mn matrix 

[H^ : FH^ : fV : (IV-68) 

has rank n. 

In equation (IV-68) 

n = order of the system (6) 

H = observation matrix defined in this section 
F = state transition matrix defined in Chapter II 

Performing the indicated multiplication for the H matrix used 
during the burn will yield a 6 x 36 matrix. To see if its rank is the 
necessary 6, it is postmulti plied by its transpose and the determinant of 
the resulting 6x6 matrix is found. For the system to be observable, 
this determinant must be non-zero. The procedure for the coast phase is 
essentially the same with the exception that application of (IV-68) will 
yield a 6 X 18 matrix. The results of this test showed that, for both 
phases, the determinant is non-zero and thus the system is observable. 

Simply because the system is observable at one point in the flight 
does not mean that it will always be so. This is due to the time-varying 
nature of both the F and H matrices. Therefore, the observability 
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condition must be checked at each state update. Thus far, no problem 
with a non-observable system has been encountered at any stage of the 
flight. 

D. Propagation Time Considerations 

The effect of aT, the propagation time, on the measurement system 
will be the final topic considered in this section. The value of AT 
during the coast is not a very critical factor because of the slowly 
varying nature of the measurements. The effect during the burn, however, 
while arising due to a rather subtle reason, is much more marked. 

The reason Ties in the measurement devices themselves. The 
accelerometers measure on-board accelerations only, those provided by the 
engine. The acceleration due to gravity, which is a function only of 
position, must be added. As will be recalled from the chapter dealing 
with the vehicle dynamics, the accuracy of the value obtained from the 
gravity model is a direct function of the propagation time. Thus the 
proper choice of aT is necessary to achieve accurate measurement results. 

E . Other Measurement Techniques 

The burn and coast measurement combinations are not the only ones 
possible or even contemplated. Some of the other possibilities include 
electron beam imagers [19], landmark tracking [20], tracking and data 
relay satellite system (TDRSS) [21] among others. Each requires accuracy, 
cost, and reliability trade-offs as well as depending on the type of 
mission (low earth, geosynch, etc.). The accelerometers-only during burn 
and horizon sensor-sun seeker during coast are reasonable combinations, 
however. 




97 


V. FLIGHT SIMULATION ANALYSIS AND RESULTS 

In this chapter actual simulation results and an analysis of them 
are supplied in order to yield an indication of what type of overall 
response- tracking characteristics could be expected to be encountered. 

In most cases data is supplied in terms of lumped position and velocity 
results with one axis of either x and x, y and y, or z and z also included 
It should be kept in mind that even though one axis (x, y, or z) may 
exhibit favorable response, the oth^r axes may not and this is the purpose 
of the lumped position and velocity tracking results. 

A. Basic Mission Description 

The primary mission for the Space Tug is given in Chapter I. This 
mission profile was chosen as a base mission for the project. In the 
studies which follow in this chapter the first burn and first coast 
phases of the Space Tug's maneuver from low earth to geosynchronous orbit 
are investigated. It was decided that these two phases of flight would 
cause the most difficulty for the estimation scheme and other burns and 
coasts would be an extension of the first. 

B. Filter Checkout 

The logical first stage in checking the entire simulation program was 
a verification of the operation of the filter itself. This checkout was 
accomplished in three steps in order to verify that the simulator was 
functioning properly. 
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The fi'rst of these tests was a deterministic one in which actual 
(perfect) state information was made available as measurement data. In 
addition to the perfect measurements, the initial estimated states were 
made equal to the actual ones. 

Next, offsets were added to the initial states but the measurements 
were kept perfect. This was to take into account the practical problem 
of initial state indeterminacy. Finally, noise was added along with 
the initial offsets to simulate the combined effect of system noise and 
measurement errors. The noise was computed during each cycle by using 
a pseudo-random number generator with mean of zero and standard deviation 
(a) a percentage of the actual states. By using the actual states as 
measurements, the problem of initial error propagation discussed in 
Chapter IV in association with indirect measurements was avoided. This 
is not the method to be used in the final program, but suffices to check 
the filter by assessing its sensitivity to noise. The variable standard 
deviation allows some parametric data to be gathered regarding the 
amount of noise tolerable in the system in order to meet design specifica- 
tions. 

The results of these tests are shown in Figure V-1 through V-3. Each 
figure contains plots to show the results of each test on the normed 
position and velocity errors. These error norms were calculated by using 
the equations below 


Xe(t) = ^(X-X)^ + (Y-Y)^ + (Z-Z)^ 
Xe(t) = ^(^X)^ + (Y~Y)^ + (Z-Z)^ 


(V-1) 

(V-2) 
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where Xe and Xe are the normed position and velocity errors respectively 

• • • 

X, Y, Z, X, Y, Z are the actual states 

/» /V /V /V /» yv 

X, Y, Z, X, Y, Z are the estimated states 


In addition, each run was performed with a propagation time (At) of 
two seconds, and P, Q, and R matrices defined as follows: 

10® 0 0 0 0 0 

0 10 ® 0 0 0 0 

p , ;o 0 10® 0 0 0 ,„.3) 

0 0 0 10° 0 0 

0 0 0 0 10® 0 

0 0 0 0 0 10 ° 


r 2500 0 0 0 0 0 ' 

0 2500 0 0 0 0 

0 0 2500 0 0 0 

0 0 0 2500 0 0 

0 0 0 0 2500 0 

0 0 0 0 0 2 500 

10 0 0 0 0 0 

0 10 0 0 0 0 

0 0 10 0 0 0 

0 0 0 10 0 0 

0 0 0 0 10 0 

0 0 0 0 0 10 


(V-4) 


(V-5) 


The results for the deterministic run showed that both the normed 
position and velocity errors remained very nearly zero throughout the 
first burn. No significant deviation was encountered. 
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Adding an initial error to the estimated states resulted in Figure V-1. 
These tests show that the filter converged regardless of the magnitude 
of the offsets. The offset had been filtered out after one, or at 
the most two, cycles. 

The results of adding noise to the measurements and errors to the 
initial state estimate are shown in Figures V-2 and V-3. The various 
curves represent different values of the standard deviation of the noise 
ranging from .01% of the actual states in (a) to .3% in (d). Once 
again, the amount of offset was not a critical factor. From these 
curves it appears that a noise of up to .1% can be tolerated and still 
fulfill the design requirements. 



Nortned Position Error (km) Normed Velocity Error (m/sec) 
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KEY 

INITIAL OFFSET (X, 1 ) 

35 35 35 

40 20 40 

20 10 20 
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Figure V-1 . NORMED POSITION AND VELOCITY ERRORS 

FOR RUN WITH INITIAL OFFSETS AND NO NOISE. 



Normed Position Error (m) 



Normed Velocity Error (m/sec) 
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Time (sec) 


Figure V-3. NORMED VELOCITY ERROR FOR 
RUN WITH INITIAL OFFSETS 
AND NO NOISE. 
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C. Accuracy of $(k+1,k) Propagation 

Figure V-4 shows position and velocity error norms during the first 
burn phase. These error norms are based on the difference between the 
actual and the estimated state vector and calculated by (V-1) and (V-2), 

The actual state vector, is obtained by direct integration of the 
state vector at the previous time instant while the estimated state 
vector, X, is obtained by inputing the present actual state information 
into the filter to get an estimated update of the state vector. This 
tells the accuracy of the one step propagation. 

Since the actual information is provided as perfect data, i.e. no 
initial offsets or noise, this comparison between the two sets of state 
information allows an indication of how well the filter is operating. 

This mode of operation is termed as the "passenger" mode since it carries 
the filter only as a passenger for each one step propagation. 

As seen from the figure, both error norms decrease almost logarithmically 
during the first half of the =240 second burn period. During the final 
half of the burn phase both show a moderate increase. However, during the 
entire burn they both stay well within tolerable limits. 

These results show that the dynamic model incorporated in the filter 
is an acceptable representation of the actual system's dynamics. However, 
since the filtered estimates do not actually drive the filter in this 
mode of operation, these errors norms give no indication about overall 
simulator package accuracy. 


Velocity Error Norm (xlO m/sec) Position Error Norm (x 10 
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Figure V-4. POSITION AND VELOCITY ERROR NORMS 
OF ONE-STEP PROPAGATION. 
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D. Covariance Init i alization 

In theory, the diagonal elements of the covariance matrix P are 
directly related to the uncertainty in the filter's output. The question 
of covariance initialization arises from the non-linear nature of the 
system. Kalman filter theory states that, for a linear system, the 
elements of the covariance matrix will get smaller as the estimated 
values approach the actual ones. However, the theory does not always hold 
for non-linear systems. The problem here is that the elements of P 
may get unjustifiably small meaning that they no longer reflect the 
error in the state vector. This causes the filter to use the one-step 
propagation at the expense of the measurements to arrive at the 
estimated state, In other words, the filter believes it is closer to 
the actual state than is actually the case. The amount of correction 
is therefore small and filter divergence can take place. The proper 
initial values for P must be chosen so that they reflect the initial 
uncertainty and yet be large enough to prevent divergence of the form 
just described. 

Although there is no straightforward analytical way to find these 
initial values, a starting point can be found by using the estimated 
uncertainty in the initial states. These uncertainties are approximately 
1000 meters in each of the X, Y, and Z directions and 1 meter/second 
for the velocity components. Since these values are standard deviations 
and the elements of the covariance matrix are variances, each must be 
squared. In an attempt to prevent divergence, the elements were initially 
multiplied by a factor of ten to yield the initial simulation test data. 
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Several runs were made with various initial values for the 
covariance matrix. The .tests served a dual purpose; to find proper 
initial covariance values, and to check the sensitivity of the system 
to changes in those values. Runs were made using values along the main 
diagonal ranging from 10^ to 10^^ with each position state and each 
velocity state uncertainty equal although the two groups were not 
necessarily equal to one another. 

The results showed that the system is not particularly sensitive 
to the covariance initialization. Although divergence did occur for 
some of the small values, the system response was essentially the same 
for initial values several orders of magnitude apart. As a result of 
these tests, the initial covariance matrix vtis chosen to be 


10^ 

0 

0 

0 

0 

0 


' 0 

io^ 

0 

0 

0 

0 


0 

10 

10® 

0 

0 

0 

(V-6) 
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0 

10^* 

0 
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0 

0 

0 

0 
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0 

0 

0 

0 

10® 
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E. Parametric Data Using Constructed States 

The tests run in this section were much akin to the ones run in the 
section dealing with the initial filter checkout, i.e., the filter was 
tested to see its response to various initial offsets and measurement 
noises. The difference was that here, states that were used as the 
measurement data during the burn were constructed from actual accelerometer 
data as opposed to the previous use of actual state information. The 
measurements were obtained by propagating the vehicle's acceleration to 
determine velocity and position data. The defining equations are: 


X = {X + y) At + Xq (V-7) 

A4. 2 

X = (X + y) {^) + X At + Xq (V-8) 


where 


X = vehicle's new velocity 
X = vehicle's new position 

X = vehicle's acceleration at time of measurement 
At = propagation time 
Xq = vehicle's initial velocity 
Xq = vehicle's initial position 
Y = noise 


A more detailed explanation of these equations is in order. There 
are two acceleration effects present here; an acceleration provided by 
the engines and one due to gravity. The on-board accelerometers only 
measure the first of these. The second must be accounted for by using 
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the gravity model previously described, with the next state estimate as 
its input. The X term includes both effects. Once again, the noise 
term was added to take care of both system noise and measurement errors. 

The filter was checked through the first two burn and coast 
phases, the first one of which is plotted here. To account for errors 
in the coast measurement, a noise was added to the coast measurements. 

Due to the angular nature of the coast measurements, the standard deviation 
of the noise was expressed in radians. When used in this manner, the 
simulation program is now in the form to be implemented on an acutal 
flight. 

As previously stated, the main objective of these tests was to 
determine the effect of noise on the filter. However, an important 
sidelight was to verify the predicted effect of initial state errors 
on the indirect measurement scheme used during the burn. Recall the 
discussion at the end of Chapter IV of these indirect measurements. It 
was stated there that since each new measurement depends on the previous 
one (through the initial condition factor), all measurements will 
include any initial error. Therefore, in general, the elimination of 
the error does not appear possible. It is enough of a problem to 
minimize it’s growth. 

All tests were carried out with a propagation time of two seconds 
during a burn and sixty seconds during a coast. Also the P, Q, and R 
matrices were the same as before. The plots shown are the errors in 
R, V, X, and X for the first burn and coast phases. Standard deviations 
were chosen on the basis of device specifications as given in [3, 4, 22, 
23, 24], 
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The plots for the first burn, shown in Figures V-5 and V-6 reveal two 
interesting features. Contrary to the results of actual states from 
Figure V-1, the initial offsets were the biggest factor in determining 
the response, while the noise, at least for reasonable values, had little 
effect. In fact for several cases, the plots for the same offsets and 
different noise were so close together as to be indistinguishable. This 
situation clearly demonstrates the effect of initial measurement error. 
While the growth of the error is not at all rapid and indeed was reduced 
slightly in some cases, the presence of that initial error can be felt 
throughout this phase. 

The estimation scheme during the coast fared better. Here the 
reduction of the error can be seen in Figures V-7 and V-8. Even 
relatively large noises did not deter the estimator from converging. 

This situation is attributable to the direct measurements that are 
available during the coast phase. 

In summary, these tests showed that the amount of the initial 
state error is a crucial factor during the burn but not the coast. 

These results were predicted using the defining equations (V-7) and 
(V-8). The slow growth of the error during a burn and its subsequent 
reduction during a coast leads to the conclusion that, given good 
initial estimates, the filter will perform as required on an actual 
mission. 
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Figure V-5. PLOT OF R AND \ 
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Figure V-6. PLOT OF X AND X ERROR DURING FIRST 
BURN USING CONSTRUCTED STATES. 
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Figure V-7. PLOT OF ERROR IN R AND V 
DURING FIRST COAST. 
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Figure V-8. PLOT OF ERROR IN X AND X DURING FIRST COAST 
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F. Study of Unfiltered Noisy Accelerometer Data vs. Filtered Noisy Data 

It was suspected that, due to the indirect measurements used during 
a burn, the filter might prove ineffectual. Accordingly, a seri'js of 
tests were made to compare the filtered and unfiltered responses. This 
was done by bypassing the filter and using the unfiltered state estimate 
as the system input. The state estimate was obtained by integrating the 
accelerometer measurements. All of the pertinent data used to carry out 
these tests was the same as that used previously. 

The results shown in Figure V-9 and V-10 bore out the prediction. 

The curves for various noises were plotted, but they were found to be 
very similar to one another. Therefore only tests using different off- 
sets are shown. 

The small differences between the two curves demonstrates dramatically 
the poor nature of the (relative) burn measurements. In essence these 
tests show that little "cleaning up" of the state estimate is possible 
during a burn due to the lack of true state measurements available. 

This shows that the six vector filter of positions and velocities with 
accelerometer data during the burn is of limited value. One possible 
remedy for this would be to incorporate a nine vector of positions, 
velocities, and accelerations. 
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Figure V-9. PLOT OF ERROR IN X AND X FOR FILTERED AND UNFILTERED 
RESPONSES 
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G. Propagation Time Study 

Filter accuracy depends on many parameters associated with the filter 
One of the basic parameters is the propagation time for the filter. At. 
This is the time interval in which the filter updates or propagates the 
state vector. In this study various values are used in the filter to 
determine a range of propagation times in which the filter can still 
function properly, as well as to determine a best or optimal propagation 
time so that errors in the state vector will meet the desired tracking 
specifications. From previous knowledge and literature on Kalman 
Filtering it was known that the propagation time for the filter should 
probably be between .5 and 2 seconds [25]. It was believed that the 
smaller the propagation time the better the filtered states should 
track the actual states and that for larger propagation times, the filter 
would be less accurate. These ideas seem reasonable when viewed from 
the standpoint of taking measurements. If measurements are taken at 
smaller time intervals, the filtered states should be closer to the 
actual states. 

The simulation runs made for this study were for three different 
propagation times, .5, 2 and 5 seconds, with .5 and 2 representing lower 
and upper bounds for the propagation time. The 5 second run was used to 
determine the rate at which the filter would diverge for an unreasonably 
large propagation time. In all cases the initial offset in the state 
vector was selected as aX^ = [1000 500 1000 212]. This offset was 
used since it was specified as a normal uncertainty in the state vector 
prior to the first burn. The measurement system used for all simulations 
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was that of accelerometer data being used to construct measurement 
states. The measurement noise was set at .01% of the vehicle accelera- 
tion. 

The results of these studies are presented graphically in Figure V-11 
and V-12. Figure V-11 shows curves of the position error norm and 
velocity error norm vs time for the first burn. These results are as 
expected. In all cases the rate in which the normed error increases is 
reduced for a smaller propagation time. The most significant difference 
in the three simulations appears in the normed velocity error curves. 

The normed position curves show a smaller amount of difference. This is 
due to the method of constructing the states in the measurement system. 
The filter is not able to reduce these errors below the value of the 

initial offset due to the lack of true measurements. The filter is 

forced to use the initial state vector as accurate. Therefore, if the 
filter tracks closely to the initial value then it is assumed that if 
a better initial value of the state vector were known the errors could 

be reduced. The error in the X and X states is shown in Figure V-12 

as an illustration of a typical single axis error response. Again, 
results were as expected. 
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Figure V-Tl. NORMED POSITION AND VELOCITY ERROR vs. TIME 
FOR DIFFERENT PROPAGATION TIMES 
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Figure V-12. ERROR IN X AND X vs. TIME 

FOR DIFFERENT PROPAGATION TIMES 
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H. Variation of Initial State Vector Offset 

As stated previously, the initial state vector prior to a burn 
will always contain an error. In GSP this initial error was modeled 
as an additive offset to the actual state vector. That is 

Xi(0) = )^(0) + (V-9) 

where i and a signify initial and actual, respectively. Therefore from 
(V-1 ) a positive number in the vector means the initial state 

is larger than the actual state. However, when the error between the 
actual and estimated state is computed a negative quantity would be 
found. The error in the state vector is computed by 

)L - X = iL (V-10) 

The purpose of this study was to observe the filter error response 
to variations in the initial state offset vector. It was decided to 
keep the initial normed position and velocity errors constant and only 
change the sign of certain values in the initial state offset vector. 

The magnitude in the initial offset vector were set at reasonable 
maximums that were specified. In the previous studies the initial state 
vector was greater than the actual state vector in each individual state. 
By constructing the initial state vector in this manner the normed posi- 
tion and velocity errors were never reduced below the initial value. 
However, in this study, where the initial state vector is greater in 
some states and less in others than the actual state vector, it is not 
clear whether the normed position and velocity errors will follow a 
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similar trend. Since the measurement system employed in the filter does 
not use exact state information it is expected that the error trend 
should be similar. 

In the simulations made for this study the measurement system was 

A 

that of the constructed measurement states and using ^(k+l/k) in the 
gravity computation. The measurement noise was again set to .01% of the 
vehicle's acceleration. The propagation time was set to .5 seconds. The 
results of the simulations are shown in Figures V-13 and V-14. 

From these results it can be seen that for some offset variations, 
these particular errors follow the same trends as the results from 
previous studies. However, for other offsets the normed position and 
velocity errors are actually reduced below their initial value. This can 
only be attributed to the particular offset value used and no general 
statement can be made regarding the reduction or increase in the error. 
Again, these point to the fact that the filter response, either good or 
bad, depends highly on the accuracy of the initial offset vector. 
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Figure V-13. NORMED POSITION AND VELOCITY ERROR vs. TIME 
FOR DIFFERENT INITIAL OFFSETS 
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Figure V-14. 


ERROR IN X AND X vs. TIME 
FOR DIFFERENT INITIAL OFFSETS 
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I. Gravity Computation Studies 

In the filtering process, a gravity model is used to obtain the 
acceleration of the Space Tug due to gravity. The acceleration due to 
gravity is then added to the sensed acceleration received from the 
accelerometers to give the total vehicle acceleration during the 
burn phase of flight. This is then used in the measurement system to 
obtain constructed measurement states. The purpose of the study made 
in this section centers around the computation of the gravity model . 

The gravity model needs as one of its inputs a value of the state vector. 

In previous simulations the value of the state vector used was the value 
X(k/k). It was realized, however, that at the time the gravity model 
was called for in the measurement system, the one step propagated state 
vector X(k+l/k) was available. The one step propagated state vector is a 
predicted value for the state vector at time given measurement in- 
formation at time t|^. By using ^(k+l/k) instead of ^(k/k) in the gravity 
model computation, it was believed that the acceleration due to gravity 
could be improved. This implied that improvement in the filtered state 
vector at time X(k+l/k+l), could be obtained. 

The simulation runs made to verify the above assumptions used an 
initial offset in the state vector of aX^ = [1000 500 1000 212] and a 
measurement noise of .01% of the vehicle acceleration. The results of 
the simulation runs are shown in Figures V-15 and V-16. The normed posi- 
tion and velocity curves both show a reduction in error as expected. 

The larger amount of difference again appears in the normed velocity 
error curves for the reasons discussed in the previous section. Figure V-16 
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shows the errors in the X and X states. There is a small improvement 
in the X state error for the two cases; however the improvement in the 
X state shows a reduction below the error of the initial X state. The 
increase in the normed error shows that the other states are increasing 
at a rate greater than the rate by which the X state decreases. 



NORfCD VELOCITY ERROR (m/sec) NORMED POSITION ERROR (km) 


128 




Figure V-15, NORMED POSITION AND VELOCITY ERROR VS. TIME 
FOR DIFFERENT COMPUTATION TIMES 
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Figure V-16. ERROR IN X AND X VS. TIME 

FOR DIFFERENT COMPUTATION TIMES 


VI. FINDINGS, CONCLUSIONS, AND RECOMMENDATIONS 


In this chapter are outlined the overall results of the Space Tug 
Navigation study. In addition are included recommendations based on an 
extension of the work reported herein using as a basis the computer 
simulation package developed under the contract. 

Based on the study of the simulation results presented in Chapter V 
and previous theoretical analysis the following findings are presented: 

1. Overall tracking accuracy is independent of the choice of the 
initial covariance matrix, P(0/0), for a reasonable range of 
values (See Chapter V-Section D). 

2. Position and velocity error norms are highly dependent on initial 
state indeterminacy. In general, the rate of growth of the error 
is a function of the individual state component inderminacies. 
(Xe, Yg, Zg, ’^e* ^e) (Refer to Chapter V-Section H) 

3. It appears as if the position and velocity error norms are re- 
lated to the state estimation filter propagation time. At, by 

II Xe(t) II : k At + II &«» It 


where, k is an unknown constant and Xg(0) is the initial state 
indeterminacy. 

4. The overall tracking accuracy is directly related to the accuracy 
of the state vector estimate used in the gravity model during 

the burn phase acceleration computation (See Chapter V-Section I). 

5. It was found in the computation of the partial derivatives of 
DCX, DCY, and DCZ that these partial derivatives were insensitive 
to a large number of intermediate partial derivative computations. 
From simulation studies of various mission profiles it was 
observed that the 


3K, 3K 8K- 3K« 

— and — where c,- is the i^^ element of the state 

’ ar_. Sc.. sc.: 1 


vector, could be set to zero with minimal degradation of overall 
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tracking accuracy (See IGM Guidance Law in Chapter III). Be- 
cause of the interrelationship of the aforementioned partial 
derivatives with others, 84 partial derivative computations 
were therefore eliminated, thus greatly reducing the computational 
burden. 

6. As shown in the results of Chapter V-Section F, the unfiltered 
noisy accelerometer results are only slightly worse than the 
filtered measurements during the burn. This was as expected 
since only relative state vector measurements are available 

at this time. It is felt that a ninth order filter formulation, 
including vehicle accelerations as states themselves, would 
dramatically improve the filtered results at the expense of a 
larger computational burden. 

7. It appears that the tracking limits specified can be met by 
the filtering scheme developed. However, the tracking accuracy 
is close to the upper limits specified and further work is 
needed to improve them. 

From the foregoing findings and the overall results from Chapter V, 
the following general conclusions can be stated: 

(1) The overall tracking accuracy can probably meet specifications. 
During the burn, however, the initial state indeterminacy 
principally determines the tracking error. The propagation 
time also affects accuracy. 

(2) The general simulator package GSP is a flexible, efficient 
simulator that is capable of handling a variety of nonlinear 
state estimation jobs. Due to its design flexibility, it 

is also useful with aircraft and rocket navigation problems 
in addition to that of the Space Tug. 

Many ideas suggest themselves as possible areas of future work. 
Although these ideas are quite varied and interrelated, they can be 
lumped under the following six major recommendations for future work: 

1) Investigation of trade-offs between navigational accuracy and 
computational requirements. The low-cost nature of the mission 
dictates the need for a low-cost, simple on-board computer. By 
reducing the volume of calculations, the computation time and 
hence the computational burden can correspondingly be reduced. 

To implement computational simplifications, the sensitivity of 
the overall tracking accuracy to each on-board calculation must 
be determined. Since the bulk of the computations are partial 
derivatives used in the one-step propagation, a good starting 
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point would be to try to further reduce these by elimination 
as described in Chapter III. In addition, a simplification 
may be possible in the gravity and drag models. 

2) Study the use of different measurement systems. As explained 
in Chapter IV, the indirect measurements used during the burn 
make reduction of the initial error difficult. Thus, it would 
be desirous to obtain direct state information from the 
measurements. One possible suggestion is to use the Interfero- 
meter Landmark Tracking (ILT) System during both the burn and 
the coast phases. 

3) Find a way to reduce the initial state error. This ties in 
closely with (2) in that one possible remedy might be found 
by using different measurement configurations. Another area 
to be investigated is the frequency, type, and accuracy 
requirements of ground based updates. 

4) Augment the present six-state vector to include the accelerations. 
By yielding better data on the accelerations, this nine-state 
vector might improve the position and velocity estimates. Such 

an approach would involve modeling the dynamics of the third 
derivative of position (jerk), however. 

5) Determine the relationship between the propagation time and 
tracking accuracy, the purpose being to determine the longest 
propagation time possible and still meet mission tracking 
specifications . 

6) Investigate the use of a non-linear filtering scheme to replace 
the present linearized one. Although it would increase the 
analytical complexity, it is expected that this new scheme would 
reduce the computational burden and increase tracking accuracy. 
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APPENDIX A 


$ MATRIX PARTIAL DERIVATIVE COMPUTATION 

In the development of the state estimator, it was necessary to 
obtain discrete linearized equations based on the system dynamics in 
order to propagate the state vector. 

The basic system dynamics considered for the 3-degree of freedom 
flight are 


X 


J DCX + X 
M g 

Y 

- 

J DCY + Y 
M g 

Z 


J DCZ + i 
LM g_ 


where T = thrust (fixed) 

M = vehicle mass (time-varying) 

Xg, Yg, Zg = gravity terms 

DCX, DCY, DCZ = direction cosines of engine gimbals 
The state vector is defined as 

f = [X Y Z M i] (A-2) 

In terms of the state equations, (A-1) can be written as 

% = f(^,t) (NON-LINEAR SYSTEM) (A-3) 

Expanding this in a Taylor Series gives 
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X(k + 1 ) = X(k) + i At + i ^ + H.O.T. 


(A-4) 


where At is the discrete propagation time and H.O.T. represents the higher 
order terms in the series expansion. Since t) and 


- . 

A ■ ^ A.» (A-4) can be written as 

9f 2 

X(k+1 ) = X(k) + f (X,t)At + ^ f ^ + n(l<) 


(A-5) 


where _n(k) is lumped noise due to series truncation. 
After defining the terms 


JACOBIAN MATRIX = A(X,t) = 


ix 


X(k) 


(A-6) 


and 


$(k + 1, k) = [A(X.t) — + I] 


(A-7) 


equation (A-5) can be written as 


X(k+1) = X(k) + $(k+l, k)f{j^,t)At + n(k) 


(A-8) 


As evidenced from (A-6) and (A-7) $(k+l, k) is found from taking 
successive partial derivatives of the system dynamics with respect to 
the elements of the state vector. It was assumed that OCX, DCY, DCZ, 
and the gravity terms are the only terms in (A-1 ) that are functions of 
the states, therefore the partial derivatives of these first three 
terms are needed in terms of known quantities. This enables propagation 
of OCX, DCY, and DCZ from one discrete time instance to the next. However, 
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generation of these partial derivatives proved to be quite laborious 
due to the number of terms involved in the I.G.M. law. 

The complete derivation of these terms with respect to one element, 
X, follows. This derivation is shown in successive steps in order to 
show the precise procedure and all quantities are as defined in I.G.M. 




sin Xp cos 
sin x; 

cos X* cos 



K = [$j]BA , BA 


^11 ^12 ^13 
^21 ^22 ^23 
^31 ^32 *^33 
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^21 

C 22 

^23 


^31 

^32 

^33 

OCX = 

sin Xp cos 

x' + 

DCY = C 21 

sin Xp cos 

X * 

DCZ = C3^ 

sin Xp cos 

X * 


nz 13 — '^y 

-22 Xy + C 23 cos Xp cos Xy 

^32 ^ ^33 ^p ^y 

Cn ■•• ^23 are functions of the states; however, through simulation 
they were shown to have little effect in the following partial derivatives. 
Hence, their derivative w.r.t. the states are considered zero. 


STEP 1. 


K . . . 

8 X ''11 ^y ^p lx ■ ^p ^y 3X 
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" ^11 
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9X 
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cos Xy sin Xp 


+ C^2 Xy - C^3 
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cos x; sin Xy -^ + 


STEP 2. 


9x' -9K^ 9K2 9 x 

P = L 4. — £ xf + . —cL 

3 X 9X 9X ax 


9x; - 9 K 3 9K. 9x 

9X 9X 9X * 9X 

9x' 3x‘ 

Now solve for — ^ and in known quantities. 
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STEP 3. 



STEP 4 . 

= lx ((«VNP)(« - «T " ^ %CgT/ - [MS-NQ!sin Xp)) 

= [c - Cj + iTg + %5gTg2 - MS sin Xp + NQ sin Xp] |)^MBy-NP) 
" [«By-NP] j|y 5 . ^ + |5j.(5Tg) * |x(’5gTg2) 

- lx (MS sin Xp) + lx (NQ sin ^p)j 
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h <Vp) = ''p lx <“p) ^ ®p lx <%) 
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lx isTg) “ j lx (Tg) - lx (5) 
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^ *p *^3 lx*"'" Xy’ * 


(sin x„) It (Ko) 


aX '‘'3' 

a 


|y(NQ Sin Xp) = N jq |y(s 1 n Xp) + sin Xp |x«!)j 


) 


+ (Q sin Xp) ^(sin x ) + sin x ^(K4) 


|x(MBj,) = H jif i- |x<TpV^^L)j + ji;(cos x^) + K, |y(sin Xy) 


-a 


■* Xj,) Sf <S> 


a 

ax 


<"P' “ " j-y^ex Ix'^) " |x(Vpx)l - |x<J> * J lx<^» 


+ P K4 lx (sin xy) * sin Xy ^ (K4) 


lx<"\> = " i''ex IxC-) " L |x(V^x> i ^ 


V fx(“V><y> * S fx<^" *y> * <='" *y l)(<''3» 


Ix(NBy) = N 


+ ^{j V L) 

ax aX'‘g''ex‘-^ 


+ By K4|x(sin xy) + sin Xy 
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Ix(ms-nq) 


= (" lx<S> ^ S |j(M)) - (n |y(Q) * Q |j(N)) 


I? ' ><y> h V> ><y) lx '"s' 

ax' ' ^^4 W IsT ^*“4* 


lx ><p> ' lx '='>* 


l"""(f)l 


= -Sin 


jtan"^ 



(Ae)fx(Al) - (A5) IxCa?) 
2 2 
(a|) + id) 


ax L ^'^21 X + l<22 Y + k23 Z] 


*^21 ^ ax ^*^21^ ^ aX ^*^22^ ^ aX ^“^23^ 


lx TJ = n It (TJ + T„ .^T (n) 


aX ''g' 'g aX 


lx 


ax ^'^21 ^ *^22 ^ ^ hs 


0 0 0 


Ix^^ (V^^^g'lx S) 



0 
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0 
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V T SL + 

''ex g 3X 


C/ \ sT 9 T bT 

\(-^\ . L — 9 . 

|\t-T y ax ax ax 


i \) 


lx = 


H) 


— = V 
ax ''ex 


! (-0 


aT aT 

._g. g. 

aX aX 


= Mi •; 9 at 
aX aX ■ aX 


9AT , ^ 

g "lx ■ 


M ..Mi •■ at 

ax ax “ ^g ax " ~§x 


ax “ aX ^ " "'g lx lx ^’^g^ 


a r -2 ^ - AC 11^ + A? M^ 

^[AC + AC ] = M 9X_ 

. .2 . 2 > 


[AC + AC ] 
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STEP 6. 


Xy> 


= COS X. 


1 


1 

. (AC + A? ) 


/ _ L-ue . 

\ AC^ + / \ AC^ + 


— (T = 2T — ^ 
3X ^'g > ^'g 3X 


3X 


(t) = 0 


8 

3X 


3(j)j 3<j>j 

(k-|i) = -q-,^ sin((|)y) -^ + "aX 


3(f)j ‘’H'J 

|Y{k-j2) = -q-|2 sin (<f)j) -^ + q32 cos(<|)j) 


3<i>n 


3(|)T ovj 

' ■‘'l3 “3? ‘*33 ~W 


3+n 


lx <'’t> = ° 


lx <‘=“= ®T> 


9<J)j Soj^ 

"lx 


lx(V''> 
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3<J), 


^(?) = X(-q^^ sin((j.^) -^+ q 3 ^ cos(<j>y) ) + 


34) 


uy«p d 4^X 

?(q,2 sin(*^) -^ + ,32 cos(+^)-^) + 


. °'*'T 

Z(-q^2 sin((j.^) -^ + q^3 eos(<{)y) ) 


3(j)., 




3G, 




9. 

3X 

= H 1 


(cos q 

11 " 

sin (^j < 








(cos 

<}>T Q]2 

+ sin (()y 

^32^ 

_J. + 
3X 




3Z 




sin (j)^ q 23 ) 

— 9 + X 
3X ^ 

(-qil 

sin (fly 


.. 

34’t 





,-qi2Sin$T 3X 

+ q ^2 cos <j). 




3(}>j 




Zg(-C 

h3 

3X^ 

^33 

C0S((f)y) 


3X. 


3(|)-|- 
~aX ■ 

8<!>J 

T 

3 <|)t 


3<t>i 


) + 


3X 


(sin Xp) = cos 


(S 



I 

(45) fy (45) - (45) |y (45) 

5 “-7 ' 


^3X 

{Ahf + (a?)^ 


3X ° 


3X ^*^22^ ” ° 


3X ^*^23^ " ° 


3X ° 
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9X '"g 


(XJ = 


(YJ = 


3X ' g 


— (Z ) = 
9X 


y See ADDENDUM to Space Tug Report 
' Contract NAS8-29852 pages 68-72 


3X ^ “ 2Y (?) ~ - T„ 


9X '■’T' 3X 


9T_ 3?^ 

g " 'g ^ 


3X 


(AT) = 




. / T \ 9G * 
- Tg) IX I 


9X ^ " aY 


ax 3X 


!!a.; !!l 

g 3X ^g 3X 


aX ^^g^ ~ ^ 


ax 

(-sin (j.^ + cos q 3 ^) -^ + 


9Y. 


{-sin q ^2 ■*■ ^32 <t>T^ 

aZ „ 9(f-i- 

(-q ^3 Sin + q 33 cos <j>^) -^ + X (-q^^ cos 


a<j). 


9<t>i 


'’si '^1 ‘*’T ■ '’sz ^ 


• * 


1. 



3(j)j 


+ z ( 
g 

-^13 

cos (j)^ 

■ 3X 

- q 33 sin 

9Hj 

an 


T _ 

n ^ 

3X 

aX 

ax 

1 

g 

^g 9X 


aX 


aY 

3Z 

= Q„, 

—3. 

+ a 

_a + 

a —21 

^21 

3X 

^^22 

3X 

‘^ZS 3X 


3X ^32 
9(^j 


9<J>j 
''T 


*'T 3X 
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STEP 7. 


3X 


94’t ^4’pt i \i 

cos y /o 3V W 9R\ . '^T / 

ax " '' 3X j R^ ^T 9X^ 


•s yj aV, aR, ) 

~T' ■ '^T “aX^J 

(v cos Y . '^T W 

I R^ ^T) ( 


aT aT„ aFK6 ) 

ax I 


3G 


T 

a^ 


n 3 ax 


a^T 9 Yt 9Vj 

^ = V^ cos Yj -^ + sin Yj 


aG 

aX 


.* 3Ar - .* 3Aii 


‘‘|''ex'-<2«*^^2xS ^ 


*2 *2 , 
) - (Ag + An + AC 


aL 


aV 


(w It + i — ^) 
^''ex aX *- aX ^ 


) (Vex*-)' 


- %(V 


Ml 
ex aX 


aV 


+ L —^) 

*- ax ^ 


a Cj 
"ax 


3Yi 


3V, 


= -V^ s1n t cos Y^ -jy 
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K 

3X 


1= *(-q„ cos qj, sin *t “Sx' 


9<i>i 


9(j)_ oip-j- 

+ n-q,2 cos - qj2 sin -^) 


9<(>n 


3<{>j '”<'y 

+ i(-q^3 cos (|ij - q33 sin 


9lj)n 


3 Hj 

~3X 


= 0 


In = 0 

3X ^ 


STEP 8. 


8(j)£j Xp q 3 i - Zp q 


3X 




X + z 


11 




■dy^)-, XX 
9X ^ 


51 - (xX + Yi + z 2 )(v 1^) I 


2 2 
R'^V 


3V 

3X 


= 0 


3R 

9X 


2 2 2 

[x^ + + r] 
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e^sinOj 

l+ejC0S9y 



COSGy) 


9sin0j 

3 X 

(1 + e^. 


2 dcose* 

(ej 3X ^ 



COS 0j) 


[1 + e^(ej + 2 cos 

30 J 

®T -lx 

cos 0 y) 

9T„ 

where 

Cl = B 22 + 62 ^ Aa + B 2 Q 

2 

C 2 = B ^2 ■'■ ^11 ■*■ ^10 

and B^g, B^-j, B^ 2 > ^ 20 ’ ^ 21 ’ ^22 constants determined 

elsewhere and Aa is the difference between semi -major axes 
of the initial Parking Orbit and the Final orbit. 

3 30t 

BX «T = “= -al 
3 . 

w ®T ^ ~n[ 



STEP 9 . 
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89 


j S'i’j 


8X 


8X 


« * >i(Tg + FK6 Tg) 


1 (.sin v) a - fos(y) • V 3 R.lL , !lr 
R ^ 8X '2 8X Rx 8X 


cos Yj 3 Vj Vj cos Yx 9 R' 


un l\y 

Yj 


8X 


8X 


V COS Y ^ „„„ 

—R R- 


i ci-c i T 3 FR 6 

^ + FK6 -ji + Tg -55P 


All partial derivative terms appearing in the above expression have 

8T 

been previously determined. It was noted that the term is common in 

80 X 


the derivative term expressions. Thus, if -ry can be expressed as a 
8T 8T 

function of and can be expressed in closed form, then all of the 
partial derivatives could be evaluated. After proper substitutions. 


80 J 

“8X 


can be shown to be. 


=>«T _ El + AS ^ 

8X 1 + A2 + A^ + A4 


where. 


X Q( 3 ,l) - Zp Q(l,l) 

El = -t— 5 W 2 


,RVX 


(XX + Y? + It) ^ 
(RV)^ 
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A2 = W5 


A3 = W6 


Sj cos 0 j 


T + 2ey cos 0 j + Oj 


6y sin 0 j 


(1 + + 2ej cos 0 j) 


A4 = W7 


Py sin 0 j 

(1 + Gy COS 0y)‘ 


and 


A5 = W8 + W9[2(ClFK6)Tg + C2FK6] 


Wi = V COS Y ^ ^ n 
R Ry 


W2 = yig + FK6(Tg)] ^sin Y 


W4 = Js[T + FK6(T )] -Vcos y 

g g 

Vy 

W5 = + FK6(Tg)] ^ sin Yy 


cos Yi 

W6 = i,[Tg + FK6(Tg)] 


W7 = Vy(W6) 

W8 = (%)(W1)(1 + FK6) 
W9=(%)(Wl)Tg 
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As a final task, must be expressed in closed form. 

At this point, it would be advantageous to discuss as it appears 
in the I.G.M. law. The equation for given in I.G.M. is. 


where it appears that this is merely an update of T^. But, acutally, 
Tg has been updated prior to this step and aT is actually a correction 
factor for the error in the update. 

The discrete equations are 


actual theor. computed 

Tg(k) = Tg(k) + AT(k) 


theor. theor. computed 

Tg(k + 1) = Tg(k) +AT(k) + 6t(k) 


and taking the first partial derivative, 

0 

actual theof. computed 

9Tq(k) _ 9T^k) 

9X M 9X 

Solving for the in terms of Tg 


9AT 

9X 


(^)s- te) 


!la. 

9X 


(A-9) 


From Step 7, 
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L 


* 3a|* o.;* 3An* 

— V— + ^An -Tv— 


{''ex L (2^ 


\ . .*^ .*^ .*2 

j - (ac + An + A? ) 


sV 

i' I V — + L 

2 ® wex 3X 9X 


) 


Substituting back into this equation using the partial derivatives 


« aT 

obtained in terms of T^, and solving for one obtains, 


3AT 

3X 






where. 


Afi = / A5(S1 - S2 - S3 - S4 - S5 - S6 + $7 + $8 + $9) \ ' , 
^ V 1 + A2 + A3 + A4 y ■ 


SI = M3 


/ 2 



+ 6j COS 0y 

ll + 

_ 

2e, cos 0x + 6 t 



h 6y sin 6j 


2 

[1 + + 26j cos 0^] ) 


I /2y \ Pj 

I [1 f cos 6^]^ 


54 = M6[q2i cos (j>^ - q^^ sin 4.^] 

55 = M7[q22 cos ())j - sin 4.^] 

56 = M8[q23 cos <})j - q^^ sin cpj] 
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57 = N3[q^-j cos (|)y + q^-j sin (j)^] 

5 8 = N4[q^2 'J’j I 32 

59 = N 5 [q ^2 <{)j + qg 3 sin (|)j] 





3X 

3Y 

3Z 

-9- + 

k « ■ -9 + 

ko 

9X 

*^22 3X 

•^23 3X 


3X 

3Y 

3Z 

— 9- + 

k -9- + 

k — 9- 

3X 

•^32 3X 

•^33 3X 
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M7 = 


M8 = 


M9 


N1 


N2 = 


N3 


N4 = 


N5 


(tf) 



A +^z V 

VexL / 

V 2 V 


( * + i X ^ 

Vex*- / 

V 2 g / 

f2An*\ 


IV • .1 
\ ex , 

1 2 ’ 9 / 

te) 


Vex'-/ 

x« 

+ 

x: 


/ Y + ^ Y \ 

) 

V ' 7 

(S) 

(' ^ ^ ^ 9 ) 


E2 = 


I e1(S 1 - S2 - $3 - S4 - S5 - S6 + S7 + 
I ’ 1 + A2 + A3 + A4 


and El, A2, A3, A4, and A5 are as defined previously. 
From (A-9) the final step in the derivation is 


S8 + S9) 
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Similar expressions can be derived for the partials with respect 
to the five remaining elements in the state vector. Due to the extreme 
size of the listing it is omitted from this report, however, it is 
available upon request as an Addendum. 



APPENDIX B 


DATA INPUT/OUTPUT 

In this Appendix, the form of data input needed for computer 
cards is shown and an example output given and discussed as to 
individual term meanings so there can be no confusion as to what the 
output yields. It should be noted that in all cases, double precision 
notation is required, although not all signifigant digits need be or 
can always be input or output. 

A. Input From Data Cards 

For different data, different formats are used. A listing of the 
groups of variables will be given, along with the format code called 
for. 

1. Earth Orbit Data 
ECC, ECCAML, ANM, TREF, AMAJOR 

Format is 5D13.8 - up to 8 decimals per number, 13 spaces total 
(including sign and exponent), D format, up to 5 numbers on a card. 

An example card field is shown below. 

COL 1 14 27 

.O1672726OD0O 8.30127160D-3 1 . 99083972D-7 

O.OOOOOOOODOO 1.49467000D11 

COL 40 53 
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2. IGM Initialization Integers 

ISW, NCYL, NFK6 

Format is 1615 - up to 16 different five-digit integers (only 3 
are used here) which are right-justified. An example card field is 


shown below 


COL 5 

10 15 

2 

2 0 


3. Tuq Characteristic Oata 


THRUST, SISP, WTSTAG, TF, TIME 

Format is 6013. 8 - up to 8 decimals per number, 13 spaces total, 
D format, up to 6 numbers on a card. An example is shown below 

COL 1 14 27 

1 .50000000004 4.44000000D02 5.85000000D04 

1 . 00000000001 4 . 1 9755900004 

COL 40 53 


4. Exact Tug State Vector 

X, OX 

Format is 6013.8 - same as is subsection 3. An example is shown 

below 

COL. 1 14 27 

-8.7143370004 -2.4140460005 -6.6738760006 

7.73020900003 -1.7556700001 -1.0550210002 

COL. 40 


53 


66 



161 


5. State Offset Vector 

XOFFST 

Format is 6D13.8 - same as in subsection 3, An example card is 
shown below 

COL. 1 14 27 40 53 66 

1000. ODO 500. ODO 1000. ODO 2.0D0 l.ODO 2.0D0 

6. Burn-Coast Initialization Times 
DTC, DPRTC, DTB, DPRTB 

Format is 6D13.8 - same as in subsection 3. An example card is 
shown below 

COL. 1 14 27 

6.00000000D016 6.00000000D01 .500000000000 

.500000000000 
COL. 40 


7. Launch Coordinates and IGM Control Times 
AZ, PHIL, TPHIT, TCHIT 

Fermat is 6013.8 - same as in subsection 3. An example card is 
shown below 

COL. 1 14 27 

8.33911243001 2.86080000001 9.9999990006 

1.50000000001 


COL. 40 
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ICHECK 


8. Simulator Mode Selection Integer 


Format is II*- integer value is 1, 2, or 3 in column 1. 


9. End-of-Burn Orbital Parameters 
RT, VT, THT, P02, E2, AR6PER, DINCL, DNODE, TGB, RANT, WDROP 

Format is 6D13.8 - same as in subsection 3. There are two cards 
required for this data. 


COL. 

1 

14 

27 


6.69584200D06 

8.36041600D03 

1. 254501 OODOO 


7.85809400D06 

1 .75470300-01 

1.80689300D02 

COL. 

40 

53 

66 

COL. 

1 

14 

27 


2.82827000D01 

9.94007300D01 

2.42300500002 


1.71 0281 00D02 

O.OOOOOOOODOO 


COL. 

40 

53 



10. End-of-Coast Time 

TCOAST 

Format is D13.8 - one number with 13 total digits or less, only 
8 of which may be decimals. An example card is shown below 

COL. 1 

4.90717848D04 
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B. GSP Output 

This section discusses and displays representation data output so 
as to clarify the purpose of the data presented. 

On the next page is shown the first output sheet containing an 
explanation of the mode selection, i.e. deterministic, full filter, or 
passenger mode, along with the IGM data input from cards as discussed 
in the previous section of this Appendix. This is for two purposes: 

(1) to insure that the data has been read properly, and (2) it provides 
a concise listing of the simulation run values to aid in oroper analysis 
of the flight data at a later time. 

On the pages following the previous sheet are representative output 
data from GSP, as commanded by subroutines NPRINT and OUTPUT. The first 
block set of data represents initial conditions data existing before 
the simulation begins. Subheadings "Begin Burn", "Begin Coast", "End 
Burn", "End Coast" appear as necessary to clearly dileneate the flight 
phase the printed data is compatable with. Note that TIME, listed at 
the beginning of each block of data, represents the reference time 
at which the data is available. This time is referenced to the initial 
Earth launch of the Space Shuttle - Space Tug combination. For instance, 
TIME = 41980.0 means 41,980 seconds (or about 11 1/2 hours) have elapsed 
from launch initiation. A defining list of symbols for the terms in 
these data blocks is given in Table B-1 . 

The sections of output data not in blocks, i.e. "FILTERED STATES", 
represents output data from the state estimator navigation system. 



tChECK = 2 

IF lChECK=l ,TRt.E DAffl L5.EL STATE ESTIMATEk IS BYPASSED 
IF 1CHECK=2, TFt STATE rSrrNATfK IS USED 

IF IChLCK=3,THE STATE cSTlPATCR IS C ARfil EDALC'JG AS A PASSENGER AND TKLE DATA USED AT THEBEGIKMNG OF EACH CYCLE 



FCC= C. 167270-01 AT L= *F-3-"13C-C2 AMf'= C.l<39CPD-06 TRFF= O.C 


AMAJOR= 0.149A7D 12 


•t 


IGP INPUT DATA 

: 

ISW = 2 NCYLs 2 NFK6= 0 


TFRUST= 

'j.l50CCi.CEi: 

CE 

S1SP = 

C-.A^ACCCCCD 

C3 

hTSTAGE= 

U.585CCCC0D 

05 

TF= 

O.IOOOOQOOC 

02-. 

- 

• J 

TIFE = 

•j.A197o59‘jC 


X- 

-■>.P71^ 337Cr 

05 

Y= 

-C.2A1ACA60U 

06 

Z* 

-0.66738760C 

07 



XC = 

C.773C2C9CD 

CA 

Y0 = 

-U. 175S67CCD 

C2 

ZD= 

-0.105502100 

03 

DTC* 

0.60COOOOOD 

02 


. ^ .* 

DFRTCs 

, .3ccccc:.cf 

02 

CTD = 

w‘.2ClcuUCCC 

01 

DPRTB= 

O.20OC0CCCD 

02 

AZ = 

0.83391L24C 

02 



PHIL = 

V ,2E6CaCL2C 

C2 

TPF'IT = 

E.9SSSS99CD 

C7 

TCHIT= 

C.15C00CCGC 

C2 






BT = 

6.6695S'.2iC 

:? 

VT = 

:.e36C416CD 

C't 

THT = 

C.l25<i5ClC0 

01 

P02 = 

0.785809400 07 



E2 = 

L .175A7C3CC 

C 

ATGPcR- 

C.1EC6893CC 

C3 

INCL = 

0.282827000 

0 2 

NODE=.. 

0.99400730D 

02 


♦ 

A 

TGC= 

T .2423T.C50C 

C3 

RANGE= 

.■.171C281CD 

C3 

V*CRCP= 

c.o 






41^ 


THE ERRCR IN THE X VECTOR IMTIALLY IS ' i, 

U.5CCCCCD ;i C.25CCC.D C3 C.SCCCOOD 03 C.ICOCOOD 01 0.5000000 00 0.1000000,01. A*"- 

TIRE AT Er\C CF CCAST (IF ANY) SHLULC UE =.C . 't 



BEGIN eLR^l 

TIi'E C.'^I^CSCCOOC Oi 

X-13 L-.73237933C 06 Y-13 -u . 2 '.91 166 AC 06 Z-13 -0.66366A39D C7 DX-13 C. 795152730 CA CY-13 0.622137150 02 OZ-13 O.0279A971O 03 

R C.6681212AC 07 V C.759A7582D OA PATH C. 336656820 00 INCL 0.287048900 C2 NODE 0. 994436660 02 RAN 0.178011L9D 03 

GX -0.9aC08876C Cw C.26189763G 01 L'l'.X C. 163688750 01 hT C. 249351710 05 WTLPS 0.549726350 05 FLOW -0.15324069D 02 

GY :.31969A7'iC CO FNYK C.A633:659C. GC CCY 0.783001720 CC E 0.667223330 05 FLPS C.15C0CC00C 05 PHIT 0.18773460C 03 

GZ 0.686317010 ul F8ZR C.29A2^'255G OC CCZ G.9177A5320 01 ISP 0.444CCCOOO 03 F/R 0.267587230 01 TGO 0.138284620 03 

GUP C.835eS32AC G2 Ch I DP-C . lC.j66A 2C0 CC CFiIY -C.9S7C5762C Cl CHICY -C. 522900330-03 CHIPC 0. 833874960 02 CHIYC -0.997162200 01 

FKl -C.11374721G 00 FK2 -C . 169 196A 3C-C 2 FK3 -C . A95 25328C-C2 FK4 -0. 7C597 7C4D-04 CHIPP 0.9C062152D 01 CHIPY 0.90740R07D 01 
CHim 0.978341790 Cv CFK2) C.17316C39C CO CHK3) C.113A1439D GC CHITP 0.239633590 OI CHITY 0.879B3541D 01 
P C.7157567AD 07 £. u . 7 1 375C 3C0-0 1 ARG C.1P703A290 C3THU C.504548C0D OlA 0.71944243D 07PSEC 0.607300390 04 

RP 0.6679AB32C C7 RA 0. 77C936550 07 VP C. 799668970 C4 VA 0.69284243C 04 AP 0.162698270 03 AA 0.718790210 03 


Tlf'E = C.A2C8CCC0C 05 

FILTEREO STATES = 

G. 7323793219C ;= -C . 2 3^1 ic.6 36 ID C 6 -U .66366439320 07 C. 795152 73AOO 04 0.62213717950 02 0.82794976150 03 

STiTc ERROR = ACTUA'--ES TI^a f fcO ... 

-Cl. 50401053160-. 0.12AiA675AAt^C5 0. 1 85 1C967A90-C4 C. A1 1778723AD-C5 -0. 33548126620-05 -0.54458836590-04 

PGSITICN cPRCR NCR.'' = C .53707 199.98C-0A 
VELCCITY error NORM = 0.547172350-04 

EXACT REASURERENTS = 

0.73237333190 j6 -0 . 23 9 1 1663.4 10 ."‘6 -0 . 663664 3932C 07 C . 795 15273440 04 0.62213714600.02 0.82794970700 03 

AVAILABLE REASURERCNTS = 

v).73C379 ni'lC w6 -0.23911663610 C6 -0.66366439320 07 C . 795 1 52 73440 04 0.62213714600 02 0.82794970700 03 

eSTlRATEC REASURERENTS = 

C. 73237934460 .'6 — . . 2 3 71166 3e4D 06 -0.66366439370 07 C . 795 15263570 04 0.62214558590 02 0.8279634101C 03 

REASURERENTS RESIOUAL = AVAILABLE - ESTIRATfC 

-0.12f4732:35n-.;l 3 . 3 1 -78 79C2' C-C 3 v.’ . 4 6303234990-C 2 C . 98!J032C628D-C3 -0 .84398748340-03 -0.13703029600-01 



CCVARIANCt HATRIX^IS^ n-C . 18 6 U533D- li-C . 116 AC 1 29D-10 0 . 32155185C-CA-0 . 5A6853360-08 
-0.gP6U533i:-ll a.g96G323AD M-C. ig'iTegAgC-ll O. 18PA7A660-C7 C.320A3393D-04 
-0. 1 16A012gC- 10-0 . 1937694 ID-1 1 C.9960323A0 0 1-Q. 1 CFS 7 C 6 9D-C7-C • 125318080-08 
C.32155185C-OA 0.188474660-07-0.106570680-07 C . 9960 3168D 01 0 .62482 393D-07 
-0.54665336C-,:a C . 3 2C 4 339 3D-C4-0 . 1 253 18080-0 6 0. 624923930-07 C.99603158D 01 
0.266811440-0 7 0.392557610-08 0 .320489070-04 3.76S91317U-C7 0. 1277C947D-07 


0.266611440-07 
C.39255761C-08 
C.32C489C7D-04 
0.76991317D-07 
0. 127709470-07 
C. 996031580 01 .. , 


it 


THE STATE TRANS IT IC:'^ MATRIX IS 

0 . lecoccc’cc Cl o.o o.c c.icccccccn ci 

C.w C.ICCCCCCCO 01 c.c c.o 

c.c o.t C.ICCCCCCCO Cl c.c 

-0.944539210-. ) 0 .48 1 18226C-C4 C . 29 1 597960-0 4 C.99953321C CC 

m 7 65 3 490- 15-0. 279 528990-03 0 . 700223640-06-0 . 13490143C-02 
c. 564529250-0 4 0.966203390-09-0.270249060-03 3.646366230-02 


C.C 0.0 

C.ICCOCCCOD 01 C.O 
C.C C.ICOOOOOOO 01 

C. 458130940-02 0.281550530-02 
C.5732936CD 00-C. 367556960-03 
C. 114324280-02 0.974713800 00 




THE ROOTS or- THE STATE TRANSITION RATRIX 
REAL IMAGINARY 

0.999592376 CC -0 . 1 1 186264E-02 




C .99999237E 
0.58659629E 
C.96659629E 
C.se718184E 
C.96718184E 

CC 0. 11 186264E-02 
CO -C.1C3226C1E-01 
CO C. 103226016-01 
CC -C.10299139E-01 
CC C. IC299139E-CI 




CELCHX = 

o*2i-CCu\.wC0 01 

Gravity terms 





XCOG=, 

-0.128750-05 

-3.1390CC-C7 

-0.422230-36 

C.O 

0.0 


0.0 

YCCG= 

-C. 139060-07 

-0. 132860-05 

C.136C3D-06 

C.C 

o 

• 

o 


0.0 

2CCG = 

-C .422UiC-'.6 

C.14346C-C6 

0.261620-05 

L • 0 

o 

• 

o 


0.0 

NOT CCNTRC 

LIABLE CET = 

N¥ = 6 


NP *= 6 



lER * i 

CeSERVABLE 

GET = 

l.CCC''C CC 








THc CPSEHVA I ic;\ r'ATiUX H IS 


c.i:c:cj:cr. 

. I w • 0 

C • c 

c . c 

c.c 

0.0 

c. : 

^ .iCZCZ'll^C Cl 


0.0 

C . L 

c.o 



C.ICCCCCCCO 

L 1 C • 0 

0.0 

0.0 


• V, 

: .c 

C. ICCCCCCZC 

Cl c.c 

c.o 

Kj • <J 

0 • 


C*G 

C.ICCCCCCCO 

Cl c.o 

c.; 


C • c 

c.o 

c.c 

O.ICOCOCOOO 01 
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Definitions follow: 


Name 


Definition 


Filtered States 

State Error 

Position Error Norm 
Velocity Error Norm 
Exact Measurements 


State vector information actually 
available to guidance law (is in 
error) , 

e = X(k) - X(k/k) 
actual 



actual measurements without noise 


Available Measurements 
Estimated Measurements 

Measurement Residual 

Covariance Matrix 
State Transition Matrix 
DELCHK 

Gravity Terms 

XDDG 

YDDG 

ZDDG 

Not Controllable - Controllable 
Not Observable - Observable 

DET 

H matrix 


actual measurements plus noise 

results which come from filter one 
step propagation, 

Z(k+l/k) = H(k+l)X(k+l/k) 

(available-estimated) measurement 
indicates how well filter is working. 

State covariance matrix 

$(k+l,k) entries 

Integration step size 

partial derivatives of gravity terms; 
going across page, 

... 

9x ’ 9y ’ ■ ■ ■ ’ 15 ^ ’ Vf 

Headings from CONOBS referring to 
observabi 1 ity/control labil ity 

determinant of a square matrix used 
in CONOBS analysis 

Elements of H matrix; 6x6 during burn 
3x6 during coast 
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Table B-1 . 

LIST OF SYMBOLS IN OUTPUT DATA BLOCK 

Symbol 


Definition 

X-13 


Exact X coordinate of the Tug in 
Apollo 13 ECLC. (meter) 

Y-13 


Exact Y coordinate of the Tug in 
Apollo 13 ECLC. (meter) 

Z-13 


Exact Z coordinate of the Tug in 
Apollo 13 ECLC. (meter) 

DX-13 


Exact X coordinate of the Tug in 
Apollo 13 ECLC. (m/s) 

DY-13 


Exact Y coordinate of the Tug in 
Apollo 13 ECLC. (m/s) 

DZ-13 


Exact Z coordinate of the Tug in 
Apollo 13 ECLC. (m/s) 

R 


Instantaneous radius of the Tug 
from earth (meter) 

V 


Instantaneous velocity of the Tug 
(m/s) 

PATH 


Flight path angle y in (in 

degrees) 

INCL 


Present inclination angle of orbit 
(in degrees) 

NODE 


Present flight orbit node (in degrees) 

RAN 


Range angle (degrees) 

6X 


Instantaneous gravity in the X coordi- 
nate (m/s2) 

GY 


Instantaneous gravity in the Y coordi- 
nate (m/s2) 

FMXM 


Acceleration due to engine thruster in 
the X coordinate (m/s^) 

FMYM 


Acceleration due to engine thruster in 
the Y coordinate (m/s^) 

GZ 


Instantaneous gravity in the Z coordinate 


(m/s2) 



Symbol 


Definition 


FMZM 

DDX 

DDZ 

WT 

WTLBS 

FLOW 

F 

FLBS 

PHIT 

ISP 

F/M 

TGO 
CHIP 
CHI DP 

CHIY 

CHIDY 

CHIPC 

CHIYC 

P 

DDY 


Acceleration due to engine thruster 
in the Z coordinate (m/s^) 

Instantaneous acceleration of the 
Tug in the X coordinate (m/s^) 

Instantaneous acceleration of the 
Tug in the Z coordinate. (m/s2) 

Vehicle mass (kilograms) 

Weight of vehicle (pounds referenced 
to sea level ) 

Mass flow rate (kg/sec) 

Vehicle Thrust (in Newtons) 

Vehicle Thrust (in pounds referenced 
to sea level ) 

- Terminal Range angle 

Specific impulse (in seconds) 

Instantaneous Tug acceleration 
due to thrusting (m/s^) 

Time-to-go (seconds) 

Xp - Pitch burn angle (degrees) 

X - Pitch burn angle rate (degrees/ 

P seconds) 

Xy - Yaw burn angle (degrees) 

X - Yaw burn angle rate (degrees/ 

^ seconds) 

One stage computation of Xp (degrees) 

One stage computation of Xy (degrees) 

( Semi-latus rectum of instantaneous 

orbit (meters) 

Instantaneous acceleration of the Tug 
in the Y coordinate m/s^ / 
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Definition 

Instantaneous eccentricity of the 
present orbit (dimensionless) 

Argument of the perigee (degrees) 

True anomaly (degrees) 

Orbital semi -major axis (meters) 

Orbital Period (seconds) 

Perigee of instantaneous orbit 
(meters) (low point) 

Apogee of instantaneous orbit (high 
point) (meters) 

Velocity at perigee (meter/sec) 

Orbit perigee (nautical miles) 

1 n.m. = 6080.27 ft. 

Orbit apogee (nautical miles) 

Velocity at apogee (meter/sec) 



i 


APPENDIX C 


SPACE TUG NAVIGATION 
SPECIFICATIONS 


The state vector errors allowable in the initial state vector prior 
to a burn, measured in a local vertical coordinate system are: 

X - position 1 .0 Km 

Y - position 0.1 lOn 

Z - position 1 .0 Km 

X - velocity 2.0 m/sec 

Y - velocity 0.5 m/sec 

Z - velocity 2.0 m/sec 

The errors in the platform attitude angles should initially be 
less than: 


3.0 min 

3.0 

3.0 

The navigation system position and velocity errors during powered 
flight should be maintained less than: 

(excluding initialization errors) 

Position (all axes) 2 Km 

Velocity (ajl axes) 5 m/sec 

Guidance placement accuracy is bound by the following limitations: 




Position 


Velocity 

Geosynchronous Orbit 50 Km 10 m/s 

Low Earth Orbit 10 Km 10 m/s 

The operation of the navigation system should not place any 
constraints on the orientation of the vehicle. 



APPENDIX D 


SYMBOLS AND DEFINITIONS OF COMPUTER VARIABLES IN THE GSP LISTING 

On the following pages are a listing of the computer variable names 
used in GSP and the corresponding variable names implied in the IGM 
guidance law. It is included to help in the understanding and use of 
the simulator. 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DLDX 

sL/aX 

DQDXD 

aQ/ax 

DLDY 

9L/8Y 

DQDYD 

aO/aY 

DLDZ 

9L/3Z 

DQDZD 

aQ/aZ 

DLDXD^ 

9L/9X 

DUDX 

au/ax 

DLDYD 

9L/9Y 

DUDY 

aU/aY 

DLDZD 

9 L/ 9 Z 

DUDZ 

aU/aZ 

DJDX 

9 J/ 9 X 

DUDXD 

aU/aX 

DJDY 

9J/9Y 

DUDYD 

aU/aY 

DJDZ 

90/ az 

DUDZD 

aU/aZ 

DJDXD 

aJ/aX 

DAYDX 

aAy/aX 

DJDYD 

aO/a'? 

DAYDY 

eAy/aY 

DJDZD 

aJ/aZ 

DAYDZ 

aAy/aZ 

DPDX 

aP/aX 

DAYDXD 

aAy/aX 

DPDY 

aP/aY 

DAYDYD 

aAy/aY 

DPDZ 

aP/aZ 

DAYDZD 

aAy/aZ 

DPDXD 

aP/aX 

DBYDX 

aBy/ax 

DPDYD 

aP/aY 

DBYDY 

aBy/aY 

DPDZD 

aP/aZ 

DBYDZ 

aBy/ aZ 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DSDX 

3S/9X 

DBYDXD 

aBy/SX 

DSDY 

SS/aY 

DBYDYD 

3By/3Y 

DSDZ 

t 

3S/3Z 

DBYDZD 

aBy/aZ 

OSDXD 

dS/dX 

DETAX 

3n/9X 

DSDYD 

3S/3Y 

DETAY 

3n/3Y 

DSDZD 

3S/3Z 

DETAZ 

3n/9Z 

DQDX 

3Q/8X 

DETAXD 

3n/3X 

DQDY 

3Q/8Y 

DETAYD 

8n/3Y 

DQDZ 

30/ 37. 

DETAZD 

3n/3Z 

ETADX 

3n/ dX 

ZTADXD 

3C/3X 

ETADY 

3n/3Y 

ZTADYD 

3c/aY 

ETADZ 

3n/3Z 

ZTADZD 

3C/3Z 

ETADXD 

3n/3X 

XDDGX 

3Xg/3X 

ETADYD 

3n/3Y 

XDDGY 

3Xg/3V 

ETADZD 

3n/3Z 

XDDGZ 

3Xg/3Z 

DXIDX 

3?/3X 

XDDGXD 

8Xg/3X 

DXIDY, 

35 / 3 Y 

XDDGYD 

3Xg/3Y 

DXIDZ 

k/dZ 

XDDGZD 

3Xp/3Z 

DXIDXD 

35/3X 

YDDGX 

3Yg/3X 



SYMBOL 

V 

DXIDYD 

DXIDZD 

C, 

/ XIDDX 
XIDDY 
XIDDZ 
XIDDXD 

')^9DYD 

r 

XIDDZD 

ZTAX 

ZTAY 

ZTAZ 

ZTAXD 

ZTAYD 

ZTAZD 

ZTADX 

ZTADY 

ZTADZ 

XIDGX 


DEFINITION 

3C/3Y 

as/9Z 

a|/aX 

a|/aY 

al/sz 

al/aX 

a|/aY 

35/aZ 

ac/aX 

a?/aY 

a^/aZ 

ai^/aX ^ 

a?/aY 

a?/aZ 

aE/aX 

ac/aY 

a^/aZ 

aig/ax 


SYMBOL 


DEFINITION 


YDDGY 

YDDGZ 

YDDGXD 

YDDGYD 

YBDGZD 

ZDDGX 

ZDDGY 

ZDDGZ 

ZDDGXD 

ZDDGYD 

ZDDGZD 

ETDGX 

ETDGY 

ETDGZ 

ETDGXD 

ETDGYD 

ETDGZD 


aYg/aY 

aYg/aZ 

aYg/aX 

aYg/aY 

aYg/aZ 

aZg/aX 

aZg/aY 

aZg/aZ 

aZg/aX 

aZg/aV 

aZg/aZ 

arig/aX 

a-rig/aY 

arig/aZ 

arig/aX 

atig/aY 

atig/aZ 

aR/aX 


DRDXD 



SYMBOL 

DEFINITION 

SYMBOL 

\ 

DEFINITION 

XIDGY 

3Cg/3Y 

DRDYD 

aR/3Y 

XIDGZ 

sig/3Z 

DRDZD 

3R/3Z 

XIDGXD 

sig/ax 

XLRVX 


XIDGYD 

aSg/aY 

XLRVY 

3(X-X)/3Y 

RV 

XIDGZD 

a?g/aZ 

XLRVZ 


ZTDGX 

a?g/aX 

XLRVXD 

3{[^y )/a.< 

ZTDGY 

a?g/aY 

XLRVYD 

3(^/)/8Y 

ZTDGZ 

a^g/aZ 

XLRVZD 

8(|^)/si 

ZTDGXD 

aCg/ aX 

DGAMX 

3y/3X 

ZTDGYD 

a?g/aY 

DGAMY 

3y/3Y 

ZTDGZD 


DGAMZ 

3y/3Z 

ETDTX 

an-j-YaX 

DGAMXD 

3y/3X 

ETDTY 

arty/aY 

DGAMY D 

3y/3Y 

ETDTZ 

alij/aZ 

DGAMZD 

aY/3Z 

ETDTXD 

arij/ aX 

DXFDX 

3Xp/3X 

ETDTYD 

3nj/9Y 

DXFDY 

aXp/aY 

ETDTZD 

anj/aZ 

DXFDZ 

3Xp/3Z 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 




9Xp/9X 

DVDX 

9V/3X 

DXFDXD 

DVDY 

ay/3Y 

DXFDYD 

9Xp/9Y 

DVDZ 

8V/8Z 

DXFDZD 

9Xp/9Z 

DVDXD 

9V73X 

DZLXX 

9{^)/9X 


, 


Zp, 

DVDYD 

9V/9Y 

DZLXY 

9(y-)/9Y 

DVDZD 

-> 9V/9Z 

DZLXZ 

Zf. 

9(j^)/9Z 

DRDX 

9R/9X 

DZLXXD 

9(^)/9X 

^F 

Zp 

9(v^)/9Y 

DRDY 

9R/9Y 

DZLXYD 

DRDZ 

9R/9Z 

DZDXZD 

Zp • 
9(--)/9Z 



''P 

DFETX 

9<{)£t/9X 

DRTXD 

9Ry/9X 

DFETY 

8<j>Fj/8Y 

DRTYD 

9Rj/9Y 

DFETZ 

8<|)py/9Z \ 

DRTZD 

9Ry/9Z 

DFETXD 

9()>£j/9X 

DGMTX 

9Yj/9X 

DFETYD 

9(j)£j/9Y 

DGMTY 

9Yj/9Y 

DFETZD 

9(|)g-,./3Z 

DGMTZ 

9Yj/9Z 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DETX 

3Ej/8X 

DGMTXD 

avj/ax 

DETY 

3Ej/9Y 

DGMTYD 

avy/aY 

DETZ 

3E. "IZ 

DGMTZD 

3yj/9Z 

DETXD 

3Ej/9X 

DVTX 

9V-|./aX 

DETYD 

aE^/aY 

DVTY 

3V^/3Y 

DETZD 

aEj/aZ 

DVTZ 

3Vj/aZ 

DFK6X 

aFK6/aX 

DVTXD 

3V^/3X 

DFK6Y 

aFK6/aY 

D^TYD 

3V^/3Y 

DFK6Z 

aFK6/aZ 

DVTZD 

dVj/di 

DFK6XD 

aFK6/aX 

DGTX 

3G^/3X 

DFK6YD' 

aFK6/aY 

DGTY 

3 G-J./ 3 Y 

DFK6ZD 

aFK6/aZ 

DGTZ 

3G^/3Z 

DFTX 

a<{)j/aX 

DGTXD 

3G^/3X 

DFTY 

a<j)j/aY 

DGTYD 

3G^/3Y 

DFTZ 

3(j)^/aZ 

DGTZD 

3G^/3Z 

DFTXD 

3<j)^/3X 

XIDTX 

slj/sX 

DFTYD 

3<t)^/3Y 

XIDTY 

alj/aY 

DFTZD 

3(f) j/ 3Z 

XIDTZ 

3i-|-/3Z 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DRTX 

3R^/8X 

XIDTXD 

9i-,-/9X 

DRTY 

3R-p/3Y 

XIDTYD 

9|y/3Y 

DRTZ 

3Ry/ 3Z 

XIDTZD 

9i../3Z 

ZTDTX 

3?^/3X 

ZTDSXD 

3A?*/3X 

ZTDTY 

aCj/3Y 

ZTDSYD 

3A?*/9Y 

ZTDTZ 

3?j/3Z 

ZTDSZD 

3A?*/9Z 

ZTDTXD 

9^/3X 

DGDX 

3G/3X 

ZTDTYD 

3Cj/3Y 

DGDY 

3G/3Y 

ZTDTZD 


DGDZ 

3G/3Z 

XIGTX 

9Cgj/3X 

DGDXD 

3G/3X 

XIGTY 

35gj/3Y 

DGDYD 

3G/3Y 

XIGTZ 

3ig-p/3Z 

DGDZD 

3G/3Z 

XIGTXD 

95gj/9X 

DDTDX 

3AT/3X 

XIGTYD 

3ig^/3Y 

DDTDY 

3AT/3Y 

XIGTZD 

35gj/9Z 

DDTDZ 

3AT/3Z 

XIDSX 

3a|*/3X 

DDTDXD 

3AT/3X 

XIDSY 

3a|*/9Y 

DDTDYD 

3AT/3Y 

XIDSZ 

3a|'V9Z 

DDTDZD 

3AT/3Z 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

XIDSXD 

8 a|*/bX 

LXIDX 

35/ 3X 

XIDSYD 

ba|*/3Y 

LXIDY 

35 / 3Y 

XIDSZD 

3A?*/3Z 

LXIDZ 

35 / 3 Z 

ETDSX 

3An*/8X 

LXIDXD 

35 / 3 X 

ETDSY 

3An*/3Y 

LXIDYD 

35 / 3 Y 

ETDSZ 

3An*/3Z 

LXIDZD 

35/ 3Z 

ETDSXD 

3An*/3X 

LETDX 

3Ar|/3X 

ETDSYD 

3An*/3Y 

LETDY 

3An/3Y 

ETDSZD 

3An*/3Z 

LETDZ 

3An/ 3 / 

ZTDSX 

3A5*/3X 

LETDXD 

3Ari/ 3 X 

ZTDSY 

3AC*/3Y 

LETDYD 

3Ar)/3Y 

ZTDSZ 

3A?*/3Z 

lETDZO 

3Ati/9^ 

LZTDX 

3A?/3x 

DDYDXD 

3Dy/3X 

LZTDY 

3AC/3Y 

DDYDYD 

9Dy/8Y 

LZTDZ 

3A?/3Z 

noYnzn 

BOy/aZ 

LZTDXD 

3A?/3X 

DEYDX 

3Ey/3X 

LZTDYD 

3AC/3Y 

DEYDY 

3Ey/3Y 

LZTDZD 

3A?/3Z 

DEYDZ 

3Ey/9Z 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

KIPTX 

9Xp/9X 

DEYDXD 

aEy/aX 

KIPTY 

V 9Xp/9V 

DEYDYD 

aEy/aY 

KIPTZ 

3Xp/9Z 

DEYDZD 

aEy/aZ 

KIPTXD 

9Xp/8X 

DK3DX 

aKg/aX 

KIPTYD 

9Xp/3Y 

DK3DY 

aKg/aY 

KIPTZD 

3Xp/3Z 

DK3DZ 

aK^/aZ 

KIYTX 

3Xy/3X 

DK3DXD 

aK^/aX 

KIYTY 


DK3DYD 

aK3/aV 

KIYTZ 

SXy^ SZ 

DK3DZD 

a Kj/aZ 

KIYTXD 

3Xy/3X 

DK4DX 

aK /aX 
4 

KIYTYD 

3Xy/»Y 

DK4DY 

aK^/aY 

KIYTZD 

3Xy/^“ 

DK4DZ 

aK^/aZ 

DCYDX 

aCy/sX 

DK4DXD 

ayaX 

DCYDY 

aCy/^T 

DK4DYD 

ayaY 

DCYDZ 

aCy/aZ 

DK4DZD 

aK^/aZ 

DCYDXD 

aCy/aX 

DMDX 

aM/aX 

DCYDYD 

aCy/aY 

DMDY 

aM/aY 

DCYDZD 

aCy/aZ 

DMOZ 

1 

aM/aZ 



SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DDYDX 

9Dy/9X 

DMDXD 

9M/9X 

DDYDY 

SDy/SY 

DMDYD 

9M/9Y 

DDYDZ 

aDy/aZ 

DMDZD 

9M/9Z 

DNDX 

8N/3X 

DDPDXD 

9Dp/9X 

DNDY 

3N/3Y 

DDPDYD 

9Dp/9Y 

DNDZ 

3N/3Z 

DDPDZD 

3Dp/9Z 

DNDXD 

3N/3X 

DEPDX 

9Ep/9X 

DNDYD 

9N/3Y 

DEPDY 

9Ep/9Y 

DNDZD 

9N/9Z 

DEPDZ 

9Ep/9Z 

DAPDX 

3Ap/-X 

DEPDXD 

3Ep/9X 

DAPDY 

3Ap/3Y 

DEPDYD 

9Ep/9Y 

DAPDZ 

aAp/aZ 

DEPDZD 

9Ep/ sZ 

DAPDXD 

9Ap/9X 

DKIDX 

9K^/9X 

DAPDYD 

3Ap/9Y 

DKIDY 

9K^/9Y 

DAPDZD 

3Ap/3Z 

DKIDZ 

9K^/9Z 

DBPDX 

3Bp/9X 

DKIDXD 

9K^/9X 

DBPDY 

3Bp/9Y 

DKIDYD 

9K^/9Y 

DBPDZ 

9Bp/3Z 

DKIDZD 

9K^/9Z 


SYMBOL 

DEFINITION 

SYMBOL 

DEFINITION 

DBPDXD 

SBp/3X 

DK2DX 

9K2/3X 

DBPDYD 

3Bp/3Y 

DK2DY 

3K2/3Y 

DBPDZD 

3Bp/3Z 

DK2DZ 

3 X 2 / 3Z 

DCPDX 

3Cp/3X 

DK2DXD 

3K2/3X 

DCPDY 

3Cp/sY 

DK2DYD 

3K2/3Y 

DCPDZ 

3Cp/3Z 

DK2DZD 

3K2/3Z 

DCPDXD 

3Cp/3X 

KIPPX 

3x'p/3X 

DCPDYD 

3Cp/3Y 

KIPPY 

3x’p/3^ 

DCPDZD 

3Cp/3Z 

KIPPZ 

3x'p/3Z 

DDPDX 

3Dp/3X 
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This Addendum is to supplement the report on the Space Tug under 
NASA contract #NAS8-29852. The purpose is not to present a complete 
derivation of the Generalized Simulation Package, (GSP) , but to give 
insight to the approach in which the one-step propogation method was 
utilized in the IGM package to estimate the desired bum angles for the 
guidance routine. 

A complete listing of all partial derivatives and l.-jmpcd parameters 
used in the partial derivative routines of the GSP program are given. 

These equations are listed in their order of computation and should 
serve as a reference for translation of the program variable names to 
the notation of the IGM package. 

Thj.$ listing goes along with the GSP program listing, available on 
x'equest from the Department of Electrical Engineering, Auburn University. 
Together they give an overall picture of the analytical approach taken 
in the Tug navigation study. In addition, a brief discussion and derivation 
showing a step-by-step process of the partial derivative computation for 
one of the terms needed in the burn angle estimations is given in the 
Appendix of the Space Tug report. 
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